FlD-fll34  829  THE  DEVELOPMENT  OF  ft  NUMERICfiL  SOLUTION  TO  THE 

TRANSPORT  EQUATION  REPORT  1  METHODOLOGV<U)  COASTAL 
ENGINEERING  RESEARCH  CENTER  VICKSBURG  MS  R  A  SCHMALZ 
UNCLASSIFIED  SEP  83  CERC-83-2-1  F/G  8/3 


me  FILE  copy 


MISCELLANEOUS  PAPER  CERC-83-2 


THE  DEVELOPMENT 
OF  A  NUMERICAL  SOLUTION 
TO  THE  TRANSPORT  EQUATION 

Report  1 

METHODOLOGY 

by 

R.  A.  Schmalz,  Jr. 

Coastal  Engineering  Research  Center 
U.  S.  Army  Engineer  Waterways  Experiment  Station 
P.  0.  Box  631,  Vicksburg,  Miss.  39180 


September  1983 
Report  1  of  a  Series 

Approved  For  Public  Release:  Distribution  Unlimited 


Prepared  for  U.  S.  Army  Engineer  District,  Mobile 
Mobile,  Ala.  36628 


83  11  IS  054 


Destroy  this  report  when  no  longer  needed.  Do  not 
return  it  to  the  originator. 


The  findings  in  this  report  are  not  to  be  construed  as  an 
official  Department  of  the  Army  position  unless  so 
designated  by  other  authorized  documents. 


The  contents  of  this  report  are  not  to  be  used  for 
advertising,  publication,  or  promotional  purposes. 
Citation  of  trade  names  does  not  constitute  an 
official  endorsement  or  approval  of  the  use  of  such 
commercial  products. 


_ Unclassified _ 

security  classification  of  This  page  0«r«  Enrer«d> 

I  REPORT  DOCUMENTATION  PAGE 


It.  REPORT  NUMBER 


2.  GOVT  ACCESSION  NO. 


Miscellaneous  Paper  CERC-83-2 


I4.  TITLE  Cand  Subrif/*) 


DEVELOPMENT  OF  A  NUMERICAL  SOLUTION  TO  THE 
TRANSPORT  EQUATION;  Report  1;  METHODOLOGY 


READ  INSTRUCTIONS 

_ BEFORE  COMPLETING  FORM 

3  RECIPieNT'S  catalog  NUMBER 


5  type  of  report  &  PERIOD  COVERED 


Report  1  of  a  series 

6  PERFORMING  ORC.  REPORT  NUMBER 


|7.  AUTHORfa; 


B.  contract  or  grant  NUMBERTa) 


Richard  A.  Schmalz,  Jr. 

9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

U.  S.  Army  Engineer  Waterways  Experiment  Station 
Coastal  Engineering  Research  Center 
P.  0.  Box  631,  Vicksburg,  Miss.  39180 

I  1.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

U.  S.  Army  Engineer  District,  Mobile 
P.  0.  Box  2288,  Mobile,  Ala.  36628 


lO.  program  element,  project,  task 

AREA  &  WORK  UNIT  NUMBERS 


12.  REPORT  DATE 

September  1983 

13  NUMBER  OF  P  AGES 

104 


U.  MONITORIPC  agency  name  »  AOORESSC// Bl/Zaran/  from  Conirollint  Olllct)  IS.  SECURITY  CL  ASS.  fo/ lA/a  rAportJ 

Unclassified 

ISa.  DECLASSIFICATION  'downgrading 
SCHEDULE 

16.  DISTRIBUTION  STATEMENT  (el  thie  Kefiort) 


Approved  for  public  release;  distribution  unlimited. 


I  17*  DISTRIBUTION  STATEMENT  (ot  abalracf  anrarcd  In  Block  20,  It  dlttoronf  from  Boport) 


ta.  supplementary  NOTES 

Available  from  National  Technical  Information  Service,  5285  Port  Royal  Road, 
Springfield,  Va .  22161. 


19.  KEY  WORDS  (Contlnuo  on  rmvoroo  oldm  II  nocofory  and  Identity  by  block  number) 

Finite  difference  method 
Hydrodynamics 
Mathematical  models 
Mississippi  Sound 
Salinity 

20.  abstract  CCiartAMM  merererwm  »»  If  imeeeeetr  Idenlllr  By  Bloc*  nuniBar; 

^This  report  constitutes  the  first  in  a  three-report  series  and  develops 
the  methodology  for  obtaining  a  numerical  solution  to  the  transport  equation. 
The  development  of  the  vertically  integrated  hydrodynamic  and  transport  equa¬ 
tions  is  presented  in  order  to  point  out  the  assumptions  made  in  a  two- 
dimensional  approach.  A  review* of  current  finite  difference  approximation 
techniques  is  presented.  Based  upon  this  review  the  spread  time  derivative 

(Continued) 


DO  /.rW  1473  ECMTION  of  •  NOV  6S  IS  OBSOLETE 


Unclassified 

security  CLASSIFICATION  OF  THIS  PAGE  fBF»n  Pmle  Fnletrd) 


_ Unclassified _ 

SeCUWITY  CLASSIFICATION  OF  THIS  PAQgflWmi  Enfant) 

20.  ABSTRACT  (Continued). 

and  flux-corrected  transport  algorithms  were  selected  for  numerical  implemen- 
^tation  (Report  2)  and  testing  (Report  3).  The  algorithm  developed  is  to  be 
included  within  the  U.  S.  Army  Engineer  Waterways  Experiment  Station  Implicit 
Flooding  Model  (WIFM) .  The  model  will  be  used  to  predict  changes  in  salinity 
distributions  due  to  alternative  dredge  disposal  practices  in  Mississippi 
Sound  and  adjacent  areas.^ 


Unclassified 


leCoaiTV  CLASSIFICATION  OF  THIS  PAOEf»Ti«n  D»lm  Enfnd) 


PREFACE 


The  methodology  for  the  development  of  a  numerical  solution  to  the 
transport  equation  is  reported  herein.  A  numerical  solution  procedure 
will  be  developed  In  Report  2  of  this  series.  Numerical  test  results 
are  presented  in  Report  3.  The  solution  procedure  will  be  incorporated 
in  a  numerical  model  to  be  used  for  evaluating  effects  of  proposed 
dredged  material  disposal  practices  in  Mississippi  Sound  and  adjacent 
areas . 

Project  administration  and  funding  were  provided  by  the  U.  S.  Army 
Engineer  District,  Mobile.  Funds  for  publication  were  provided  by  the 
Coastal  Engineering  Research  Center  (CERC) . 

Literature  review  and  methodology  development  were  conducted  by 
the  U.  S.  Army  Engineer  Waterways  Experiment  Station  (WES)  in  the  Hy¬ 
draulics  Laboratory  (HL)  under  the  general  supervision  of  Messrs.  H.  B. 
Simmons  and  F.  H.  Herrmann,  Jr.,  Chief  and  Assistant  Chief,  respectively, 
HL,  and  in  CERC  under  the  general  supervision  of  Drs.  R.  W.  Whalin  and 
L.  E.  Link,  Chief  and  Assistant  Chief,  respectively,  CERC,  and  Mr.  C.  E. 
Chatham,  Chief,  Wave  Dynamics  Division  (WDD) .  This  report  was  prepared 
by  Dr.  R.  A.  Schmalz,  Jr.,  WDD. 

Commanders  and  Directors  of  WES  during  data  acquisition  and  analy¬ 
sis  and  the  preparation  and  publication  of  this  report  were  COL  Nelson  P. 
Conover,  CE,  and  COL  Tilford  C.  Creel,  CE.  Technical  Director  was 
Mr.  F.  R.  Brown. 


TABLE  OF  CONTENTS 


PREFACE  .  1 

LIST  OF  FIGURES .  3 

LIST  OF  TABLES .  3 

PART  I:  INTRODUCTION .  4 

PART  II:  MODEL  EQUATION  DEVELOPMENT  .  6 

1.  Hydrodynamic  Equations  .  6 

2.  Transport  Equation  .  24 

3.  Equation  of  State .  29 

4.  Compilation  of  the  Complete  Set  of  Equations  .  30 

PART  III:  LITERATURE  REVIEW .  32 

1.  Sheng  Work  at  Aeronautical  Research  Associates  of 

Princeton  (ARAP)  .  32 

2.  Fully  Multidimensional  Flux-Corrected  Transport 

Algorithms  for  Fluids  .  4l 

3.  Holly-Preissmann  Scheme  .  45 

4.  Method  of  Second  Moments .  51 

5.  Balanced  Expansion  Technique  .  58 

6.  Stone  and  Brian  Technique .  68 

7.  An  Analysis  of  the  Numerical  Solution  of  the 

Transport  Equation  .  75 

8.  The  Leendertse  Formulation  .  79 

9.  Additional  Methods  and  Considerations  .  91 

PART  IV:  NUMERICAL  METHOD  SELECTION  AND  DEVELOPMENT  .  95 

REFERENCES .  100 


2 


LIST  OF  FIGURES 


Figure  No.  _ Title _  Page 

1.  Notation  for  presentation  of  Sheng  schemes  35 

2.  Hoi ly-Preissmann  scheme  notation  46 

3.  Method  of  moments  advection  procedure 

(one  dimension)  53 

4.  Method  of  moments  advection  procedure 

(two  dimensions)  56 

5.  Chan  definition  of  variables  60 

6.  Chan  definition  for  extension  to  include  diffusion  66 

7.  Leendertse  space  staggered  grid  system  80 

8.  Development  of  a  numerical  method  for  the  transport 

equation  99 

LIST  OF  TABLES 

Table  No.  _ Description _  Page 

I  Major  Numerical  Methods  33 

II  Chan  Auxiliary  Relations  63 


3 


DEVELOPMENT  OF  A  NUMERICAL  SOLUTION  TO 
THE  TRANSPORT  EQUATION 

Report  1:  METHODOLOGY 

PART  I :  INTRODUCTION 

This  report  constitutes  the  first  report  in  a  three-report  series 
and  develops  the  methodology  for  obtaining  a  numerical  solution  to  the 
transport  equation.  The  algorithm  developed  is  to  be  included  within 
the  Waterways  Experiment  Station  Implicit  Flooding  Model  (WIFM)  [1]. 

The  development  of  both  the  vertically  integrated  hydrodynamic  and  trans¬ 
port  equations  is  presented  in  PART  II  in  order  to  point  out  the  assump¬ 
tions  made  in  a  two-dimensional  approach.  A  literature  review  of  current 
research  in  finite  difference  approximation  techniques  to  the  transport 
equation  has  been  conducted  in  order  to  determine  the  most  effective  ap¬ 
proach  for  simulating  salinity  levels  in  Mississippi  Sound.  Results  are 
detailed  in  PART  III.  The  numerical  method  selection  and  proposed  form 
of  development  is  the  subject  of  PART  IV. 

It  is  instructive  here  to  note  the  component  tasks  in  the  develop¬ 
ment  of  a  mathematical  model.  Initially,  one  must  decide  what  form  of 
the  equations  is  to  be  approximated.  Certain  simplifications  and  assump¬ 
tions  must  be  made  to  obtain  the  system  of  equations.  Secondly,  one 
must  select  a  suitable  numerical  approximation  to  the  equations.  Next 
the  numerical  approximations  must  be  computationally  implemented.  The 
efficacy  of  the  approximation  must  then  be  tested  by  comparing  simulated 
results  against  known  solutions  under  simplified  boundary  and  flow 


conditions.  Finally,  empirical  coefficients  (friction  and  dispersion) 
must  be  adjusted  during  simulation  of  measured  prototype  conditions  and 
subsequently  verified.  This  report  presents  the  results  for  the  first 
two  steps  in  this  process.  Subsequent  steps  will  be  presented  in  sepa¬ 
rate  reports. 


PART  II:  MODEL  EQUATION  DEVELOPMENT 


The  hydrodynamic  and  transport  equations  are  developed  in  two 
dimensional  vertically  integrated  form.  The  time  averaging  concepts 
necessary  in  describing  the  turbulence  are  introduced  as  necessary  in 
the  derivations  but  are  not  treated  in  detail.  An  equation  of  state  is 
developed  which  effectively  couples  the  hydrodynamics  and  transport 
equations.  Finally,  the  complete  equation  sets  are  presented  for  the 
coupled  and  uncoupled  cases. 


1 .  Hydrodynamic  Equations 

The  general  equations  of  the  classical  hydrodynamics  for  imcom- 


pressible  flow  are  given  following  Lai's  development  as  follows  [2] 


3u  8v 
3x  3y 


+ 


(1.1) 


(1.2) 


(1.3) 

(1.4) 
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u,v,w  =  Velocity  components  in  the  x,  y,  and  z  directions, 
respectively 

F^,Fy,F^  s  Body  forces  in  the  x,  y,  and  z  directions,  respectively 
p  2  Fluid  density 
p  2  Fluid  viscosity 
p  2  Pressure 
t  2  Time 

The  following  assumptions  are  made: 

1.  The  water  is  not  deep  compared  with  the  length  of  the  wave  an( 
the  shallow  water  theory  applies. 

2.  The  vertical  velocity  of  flow  is  small. 

3.  The  vertical  acceleration  of  the  fluid  particle  is  very  small 
compared  with  the  acceleration  of  gravity,  g,  and,  hence,  can 
be  neglected. 

4.  The  pressure  is  hydrostatic  (from  the  above  assumption). 

5.  The  frictional  resistance  coefficient  for  unsteady  flow  is 
the  same  as  that  for  steady  flow,  thus  can  be  approximated 
from  the  Chezy  or  Manning  equation. 

6.  Only  shear  stresses  due  to  horizontal  velocity  components  are 
significant. 

7.  The  bottom  of  the  embayment  is  rigid  or  relatively  stable  and 
fixed  with  respect  to  time. 

8.  The  water  is  nonhomogenous  but  incompressible.  The  density 
induced  flow  appears  only  in  the  pressure  gradient  terms. 

1.1  Continuity  Equation 

If  we  consider  (1.1)  and  integrate  over  the  vertical  from  the 


bottom  Zj^(x,y)  to  the  surface  r|(x>y»t)  we  obtain 


I  I 

J  9y 


dz  +  w(r|) 


w(Zj^)  =  0 


(1.1.1) 


From  Leibniz  rule  we  may  write 


'I 

/“ 


an  .  ‘'^b  . 

ax  ■  “^^b^  ^ 


au  , 
a^ 


(1.1.2) 


M  Q2  I 


Employ  the  kinematic  boundary  condition;  namely,  for  F(x,y,z,t)  =  0 
as  a  boundary  surface  assume  any  particle  on  the  surface  regains  on  it 


implies 


DF  aF  aF  ^  aF  ^  aF 


Consider  z  =  r|(x,y,t)  at  the  free  surface,  then 


F  =  z  -  ti(x,y,t)  =  0  .  Hence 


5?:  =  o=>^  +  u^  +  v^  =  w 

Dt  ^  at  “  8x  ay 


(1.1.4) 


At  the  bottom  z  = 


Hence 


Zjj(x,y)  . 


az  az, 

DF/Dt  =0=>u^+V5-^=w 

ax  ay 


(1.1.5) 
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Returning  to  (1.1.1) 


fc  /  “ 


hi 


''  "  ''«b>37 


The  sum  of  the 


The  sum  of  the 


v(n)  1^  +  w(ri)  -  =  0  (1.1.6) 


terms  is  zero  from  the  bottom  boundary  condition, 
terms  is  equal  to  dr]/dt  from  the  free  surface 


boundary  conditions.  Thus  we  obtain: 


i?"l;  J  I 


K-  I  V  dz  = 


(1.1.7) 


Denoting 


M 

/  “ 


dz  and  V  = 


one  obtains 


I?  ^  It  l(n  -  zjai  t  i-  Kn  -  zjv)  =  0  (1.1.8) 


at  ax 


Letting  h  =  rj  “  dropping  the  bar  notation  with  the  under¬ 

standing  henceforth  we  are  considering  vertically  averaged  quantities 
one  obtains 


I?  *  Ij  *  h  “ 


(1.1.9) 
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1.2  Equations  of  Motion 

2 

Consider  the  z  motion  equation  in  wi^ii-h  Dw/Dt  =  0  ,  pA  w  -  0 
from  assumptions  2,  3  and  6,  respectively.  Thus  we  obtain 

PF^.|E=0  (1.2.1) 

is  replaced  by  -g  due  to  the  following  assumptions 

1.  The  vertical  component  of  the  Coriolis  force  is  negligible 
with  respect  to  g  . 

2.  The  vertical  tide  generating  force  component  is  negligible 
with  respect  to  g  . 

If  we  integrate  the  above  relation  from  an  arbitrary  level  z  to  the 
water  surface  obtain 


p(z)  =  p(n)  + 


/ 


pg  dr 


(1.2.2) 


To  depth  average  we  consider  the  following  relations  to  hold,  in  which 
the  bar  quantities  are  depth  averages  and  the  prime  quantities  the  local 
fluctuation  from  these  averages. 


p(z)  =  p  +  p' (z)  (1.2.3) 

p(z)  =  p  +  p' (z)  (1.2.4) 


(1.2.5) 


Taking  the  partial  derivative  of  (1.2.2)  with  respect  to  x  obtain; 
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§£  _  3p(n)  9_ 

dx  dx  dx 


pg  dr 


(1.2.6) 


Integrating  (1.2.6)  over  the  vertical  we  write: 


/  = 


pg  dr  dz  (1.2.7) 


Let  us  consider  terms  ®  ■  (D  »  separately,  in  turn.  Evaluating  (T) 


invoke  Leibniz  rule 


fh 


(p  +  p' )  dz  = 


(p  +  p')  dz 


^  (5  +  p')  aT  -  (p  *  p')  1“ 


(1.2.8) 


Note  further  from  (1.2.5), 


I  I 

J  (i  +  p')  dz  =  p(n  -  zi,)] 


-  ^  P  ( 15  ■  ^ 


Thus  we  finally  obtain  for  (T)  the  following  expression: 


(1.2.9) 


j  I;  (P  *  P') 


z  )  ^  + 
dx 


p.  ^ 

P  8x 

n 


(1.2.10) 
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If  we  let  h  =  r)  -  z.  and  assume  p’  8z,/9x  ,  p'  9ri/9x  =  0  ,  then 

D  I  2^  D  r] 

we  obtain 


/  I;  ^  ^  if 


(1.2.11) 


Evaluating  note  9p(r|)/9x  is  not  a  function  of  depth;  thus  we 


obtain 


where 


Pg  =  p(n) 


Next  consider  the  iterated  integral  expression  for  Q)  as  follows: 


J  ~  ®  ■  ^\\ 


=  g  [I  (n  -  -  3  (1.2.13) 


*  If  <1  -  2)  d.  =  g  1^ 


s  ds  =  g 


(1.2.14) 


s  =  ri  -  z 
ds  =  -dz 


Observe  from  Leibniz  rule 
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gp 


•I 

I 


dz  =  gf 


8x 


fe  / 


az. 


(n  “  z)  dz 


9n 


-  (n  -  z,)  3/-  -  (n/  n) 


(1.2.15) 


Letting  h  =  (]  -  obtain: 


(1.2.16) 


and  the  evaluation  (^  is  complete.  Assembling  all  our  results  we 
obtain  finally 


h 


^  = 
ax 


h 


^Pa  - 

^-Pgh 


9£ 

^  2  ax 


(1.2.17) 


Analagously,  the  expression  for  the  y  gradient  is  given  by 


h 


^Pa  - 


2  ay 


(1.2.18) 


Thus  we  have  employed  the  z  motion  equation  to  evaluate  the  horizontal 
pressure  gradients  in  the  x  and  y  motion  equations.  The  expressions 
obtained  above  for  these  gradients  include  the  atmospheric  pressure 
anomaly,  the  water  surface  elevation  gradient,  and  the  density  gradient 
created  by  horizontal  variations  in  salinity  and  temperature. 

Let  us  next  consider  the  material  derivative  (left  hand  side)  of 
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(1.2)  the  X  motion  equation.  Expanding  the  material  derivative  and 
adding  (1.1)  we  obtain 


r  3u 

8u 

3u  . 

3u  ^ 

{  3u 

.  3v 

3w\l 

By 

"  3z 

3y 

3zjj 

^  3(u^)  3(uv)  3(uw) 

^  3t  3x  3y  3z 


3(pu)  3(pu^)  ^  3(puv)  3(puw) 

3t  3x  3y  3z 


(1.2.19) 


Integrating  the  last  result  over  the  vertical  (which  also  holds  for 
compressible  flow) 


(1.2.20) 


Again  employing  Leibniz  rule 


3_ 

3t 


/ 

2b 


(pu)  dz  = 


dz 


+  pu 


I?  - 

n 


(1.2.21) 


lA 


m 


k-v 


m 


V*.-. 

> 


>/■ 
>  ■. 


k  f 


dz 


"b 

+  puu 


^  _ 

3x 


puu 


9z^ 

9x 


(1.2.22) 


9y 


I  I 

/ 


(puv)  dz  = 


/ 


"b 

+  puv 


(puv)  dz 

9Q 

9y 


puv 


9z, 
_ t 

9y 


(1.2.23) 


n 

/ 


(puw)  dz  =  puw  I  -  puw 

n 


(1.2.24) 


Notice  in  the  above  9Z^/9t  =  0  ,  since  the  bottom  is  assumed 
rigid.  Denoting  terms  in  (1.1.4)  by  (T) ,  and  those  in  (1.1.5)  by  Q), 
we  obtain  the  following 


9_ 

9t 


n 

/ 


(pu)  dz  -  pu 


§!1  +  9_ 

9t 


(puu)  dz 


(D 


+  puu 


9x 


-  puu 


9a  +  9_ 

9x  9y 


n 

/ 


(puv)  dz 


,(D 


+  puv 


9Z, 

t 

9y 


-  puv 


% 


+  puw 
9y 


.<3)  .(D 


-  puw 


(1.2.25) 
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v  ■ 


*1 


The  above  expression  reduces  to  the  following: 


^  Cpa)  dz  t  |j  ^ 


(puu)  dz  + 


M 

-  f 

ay  J 


(puv)  dz  (1.2.26) 


The  following  terms  are  next  defined 


n  n 

hu  =  (n  -  Z^)  u  =  J'  u(z)  dzjhv  =  (n  “  Z^)  v  = 

h  h 


where 


v(z)  dz  (1.2.27) 


u(z)  =  u  f  1  +  u'  (z)] 
v(z)  =  V  Q  +  v'  (z)] 


■I  M 

=  J"  udz+  y'uu'dz  =  hu  +  u  J  u' 


(1.2.28) 


-r^ 


dz  +  /  V  v'  dz  =  hv  +  V  I  v'  dz 


(1.2.29) 


^b  ^b 


Therefore  we  must  have 


i  U-  dz  =  i  V' 


dz  =  0 


(1.2.30) 


Rewriting  (1.2.26)  and  employing  assumption  8  we  obtain 


- 


P  ^  (uh)  + 


't 

P  1^  J  u  u  [l  +  2u'(z)  +  u'(z)^]  dz 


•I 

P  1^  J"  u  V  []l  +  u' (z)  +  v’(z)  +  u'(z)v'(z)^  dz  (1.2.31) 


{l  +  [u'(z)]^}  dz  =  ^  J'  (1  +  u'(z)v'(z)]  dz 


Noting  the  depth  integral  of  the  product  of  a  bar  and  primed  quantity  is 
zero,  we  finally  obtain 


^  [It  *  lx  (3uvh)j  (1.2.32) 


Analogously,  one  obtains  the  following  expression  for  the  left  hand 
side  of  the  y-motion  equation  (1.3)  after  integration  over  depth 


p  [h.  ^  h 


(1.2.33) 


where 


■I 


-u 


1  +  [v'(z)]  }  =  r  I  [1  +  u'(z)v'(z)]  dz 


17 


Thus  P  in  (1.2.32)  and  in  (1.2.33)  is  the  same  quantity  and  is  usually 
assumed  equal  to  unity. 

Let  us  now  consider  the  right  hand  side  of  the  x  and  y  motion 

equations.  Cnsidering  the  F  and  F  terms,  we  obtain 

X  y 

F=Qv+g+G  F=-ftu+g+G  (1.2.34) 

X  X  y  '^y  y 

where  G  is  the  tide-generating  force,  Q  is  the  Coriolis  factor, 
fl  =  2w  sin  (])  (w  =  angular  velocity  of  the  earth  rotation,  <|)  is  lati¬ 
tude),  and  g  ,  g  are  the  components  of  gravity  in  the  horizontal. 

X  y 

Assume  the  following: 


1. 

8x 

’  X 

«« 

Qv  , 

2. 

,  G 

’  y 

A 

A 

A 

A 

Qu  , 

then. 


F  =  Qv  ,  F  =  -Qu 
X  ’  y 


Integrating  over  the  vertical 


'I 

/ 


Qv  dz  =  Qvh 


•I 

/- 


Qu  dz  =  -Quh 


(1.2.35) 


where  u  ,  v  are  vertically  averaged  velocity  components,  and 
h  =  n  -  Z  . 

For  a  turbulent  flow  an  eddy  viscosity  e  is  employed  in  the 


place  of  the  dynamic  viscosity  p  .  The  terms  become  using  and 

for  horizontal  and  vertical  eddy  viscosity,  respectively: 


It  is  further  assumed  the  second  and  third  terms  are  neglible.  If  simi¬ 
lar  assumptions  are  made  for  the  other  terms  in  (1.2.36)  we  obtain 


--  IS  *  s; 


(1.2.40) 


r  e  ^ 

J  "  W  ay2 


„2, 


dz  =  s, 

ax^  ay  V  ^  L  ay^ 


..  /a^v  .  a^v 
\ax^  ay^ 


(1.2.41) 


These  terms  are  retained  in  the  motion  equations  and  serve  to  stabilize 
the  numerical  approximations.  The  vertical  eddy  viscosity  terms  are  in¬ 
tegrated  over  the  vertical  as  follows. 


^sx  ^bx 


(1.2.42) 


a^v  .  /av 

e  — ^  dz  =  £  (  ^ 

az2  ^  \ 


=  t  ■ 
sy  by 


(1.2.43) 


where  t  ,  X  are  surface  stresses  and  X,  ,  X,  bottom  stresses, 
sx  ’  sy  by  ’  bx 

Consider  the  bottom  stress,  x^^  ,  as  follows: 


'b  =  <=£  0  r 

where  is  a  drag  coefficient  and  is  the  fluid  velocity.  Letting 


Chezy  c  =  J2g/C, 


obtain; 


where 

=  drag  coefficient 

=  air  density 

V  =  wind  speed 
w 

Assuming  the  shear  stress  varies  linearly  with  depth  we  obtain 


(1.2.46) 


where 

A  =  /T  +  1) 

b  s 

=  pressure  intensity  produced  by  the  wind 
s  =  distance  in  the  downwind  direction 
Integrating  (1.2.46)  over  the  vertical 


21 


n 


c;xp 
f  w 

2h 


(1.2.47) 


with 

K  =  (c^\p^)/2 

If  G  is  the  angle  between  the  wind  direction  and  the  +  x  axis, 


2 

=  KV 

cos  6 

(1.2.48) 

sx  w 

2 

=  KV 

sin  6 

(1.2.49) 

sy  w 

Assembling  our  results,  we  obtain  the  final  expression  for  the 


depth  integrated  motion  equations 


2“  2“ 

p  1^  (uh)  +  P  1^  Ou^h)  P  1^  Ouvh)  =  pQvh  +  he,  (  ^ 


9P„ 


1/2 


2"  2  "A 

P  It  ^  lx  P  1^  (Pv\)  =  -  pQuh  +  h8j^ 

\9x^  9y 


-  h  \  ^  +  pg 


+  3e  §n  +  gh  ^  )  -  pg  (u  +  v) 


1/2 


9y  2  9y 


—  +  KV^  sin  6  (1.2.51) 

w 


Setting  P  =  1  and  expanding  the  left-hand  side  of  the  above  two 
equations  one  obtains 
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3v  .  3v  ^  dv 

at  “  ax  ''  ay  “ 


-  Qu  +  e 


2  2 

a  V  a  V 

2  2 

ax^  ay 


(u  +  v) 


V  o 

+  ^  V  sin  e 
ph  w 


::ri  +  3„  ^  +  gh  9£ 
ay  ay  2  ay 


(1.2.55) 


2 .  Transport  Equation 

The  general  transport  equation  is  given  for  laminar  flow  as 


as  9s  ^  9s  .  ds  •  d 

at  “ax  ay  az  “  ^ 


(“x  t) 

{\  u)  *  h  {\  S) 


(2.1) 


where 


s  =  concentration  of  the  material  of  concern 


D  =  molecular  diffusion  coefficient  in  the  x  direction 

X 

Dy  =  molecular  diffusion  coefficient  in  the  y  direction 

D  =  molecular  diffusion  coefficient  in  the  z  direction 
z 

t,x,y,z,u,v,w  are  as  previously  defined. 

For  a  turbulent  flow,  the  eddy  dispersion  is  significantly  greater  than 
the  molecular  diffusion.  The  following  analogous  formula  holds  where 
time  averaging  over  the  time  scale  of  the  turbulence  has  been  performed. 


as  9s  9s 

—  +  u  +  v  r —  + 

at  ax  ay 


—  =  —  (k  —  ^ 
az  ax  \  X  dx  / 


{sU)  "li 


(2.2) 


where 


K^,  Ky,  and  are  turbulent  eddy  dispersion  coefficients. 
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Equation  (2.2)  may  be  written  in  conservation  form  by  arlding  s  times 
Equation  (1.1)  (namely,  zero)  to  the  left  hand  side  to  obtain 


^  a(us)  a(vs)  ^  3(ws)  -  3_  — W  —  (^K  |^W2.' 

at  ax  ay  az  ax  x  ax/  ay  \  y  ay/  az  \  z  az/ 

This  form  of  the  equation  is  then  depth  integrated.  Considering  the 


left  hand  side 


a(ws) 


(2.4) 


Each  term  is  expanded  employing  Leibniz  rule,  and  the  expressions  below 


are  obtained 


a  /  .  b  ^ 

-  /  s  dz  +  ^  s  -  ^  s 

J  /  Z 

Z,  bottom  rigid 


^  us  dz  +  (us)| 


•  \ 

J"  vs  dz  + 


^  (vs)  -  (vs) 


(2.5) 


ws  -  ws 

'  n  z, 


from  kinematic  considerations. 


The  left  hand  side  (2.4)  then  reduces  to  the  following  expression 


Let 


=  s  +  s',  v  =  v  +  v',  u  =  u  +  u',h  =  n“Z 


where 


n  n 

dz  = 


J'  s'  dz  =  f  v'  dz  =  J"  u' 

i 


dz  =  0 


Each  of  the  three  terms  in  (2.6)  is  expanded  to  yield 


at 


s 


dz  =  1^  (hs) 


9 


r 


(us) 


n 

r  (us  + 


u' s 


+ 


s'u  +  u's' ) 


dz 


where 


J'  u' s '  dz  =  h(u’ s ' 


(vs)  dz  -  J"  (vs  +  v’s  +  s'v  +  v's')  dz 


=  ly  ^(hvs)  +  h(v's')J  (2.10) 


where 


J'  v's'  dz  =  h(v's')^ 


For  the  left  hand  side  we  then  finally  obtain 


Consider  now  the  right  hand  side  of  (2.2).  As  previously  men¬ 
tioned,  the  turbulent  dispersion  coefficients  are  developed  from  time 
averaging  the  turbulence.  Specifically, 


Is  (''x  Is)  =  -  Is 


(2.12) 


^  (k  - 

9y  \  y  9y/ 


—  fv"s") 

ay  ®  ^t 


(2.13) 


^  (k  - 
az  v'^z  az  J 


—  (w"s") 

az  ®  ^t 


(2.14) 
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where  t  denotes  time  averaging  and  the  double  prime  denotes  a  time 
fluctuation.  The  terms  on  the  right  hand  side  can  now  be  integrated 
employing  Leibniz  rule  to  obtain 


n  n 


"b 

+  (u"s"). 


3Z 


9x 


n  n 


ay 

^b 

+  (v"s"), 


az 


ay 


•  I 

/ 


(w"s")^  dz  =  (w"s")j. 


-  (w’’s"), 


(2.15 


(2.16 


(2.17 


If  we  assume  all  terms  in  the  above  three  relations  are  zero  except  the 
integrals,  the  last  equation  is  removed  from  further  consideration. 
Additionally  if  we  bring  the  last  two  terms  in  (2.11)  to  the  right  hand 
side  we  obtain 

k  I  ^  ^  ^“"""^th]} 

-  L  I  h  r(v’s'),  +  (v"s").J  i  (2.18 


where 


(P  +  Pq)  (v  -  v^)  =  \  (3.1) 

where 

p  =  pressure  in  atmospheres 

p^  =  baseline  pressure  in  Jtmosphere 

V  =  specific  volume  in  ml/gm 

=  baseline  specific  volume  in  ml/gm 

X  =  constant  [(ml/gm)  atml 
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and  ^  ,  Pq  *  and  are  functions  of  temperature  and  salinity.  Eckart 

has  reported  the  following  relationships 


A  =  1779.5  +  11.25T  -  0.0745T^  -  (3.8  +  0.01T)S 

V  =0.698 
o 

p  =  5890  +  38T  -  0.375T^  +  3S 


Since  water  is  only  slightly  compressible,  we  may  neglect  pressure 
effects  in  (3.1)  by  setting  p  =  0  and  obtain 


V  A  +  V  p 
0*^0 


(3.2) 


where  p  is  the  density  in  (gm/ml)  .  In  (1.2.54)  and  (1.2.55)  p  is 
given  as  in  (3.2)  with  p  set  equal  to  l(gm/ml).  In  order  to  de¬ 
scribe  temperature  an  equation  similar  to  (2.21)  must  be  considered.  It 
is  proposed  at  this  time  to  specify  a  temperature  distribution  directly 
from  measured  data  rather  than  approximate  (2.21)  for  temperature. 

4.  Compilation  of  the  Complete  Set  of  Equations 

Relations  (1.1.9)  (Crutinuity) ,  (1.2.54)  (x-motion),  (1.2.55) 
(y-motion),  (2.21)  (Salinity  transport),  and  (3.2)  (State),  constitute 
the  complete  set  of  density  coupled  equations. 

In  cases  where  density  effects  may  be  neglected  within  the  hydro¬ 
dynamics,  the  system  assumes  an  uncoupled  form  with  the  following  for¬ 
mat.  Equation  (1.1.9)  remains  unchanged.  Equation  (3.2),  the  state 
equation,  is  no  longer  needed.  The  pressure  gradient  terms  in  (1.2.54) 
and  (1.2.55)  reduce  to  the  following  relations: 
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It  is  proposed  to  initially  consider  the  uncoupled  system  of 
equations.  In  this  manner  the  suitability  of  the  numerical  approxima¬ 
tion  to  (2.21)  can  be  assessed  directly. 
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PART  III:  LITERATURE  REVIEW 


The  transport  equation  exhibits  both  hyperbolic  and  parabolic 
characteristics.  For  convection  much  larger  than  dispersion  as  is  the 
case  in  estuarine  systems,  the  equation's  character  becomes  predomi¬ 
nantly  hyperbolic.  It  is  this  property  which  makes  numerical  approxi¬ 
mation  difficult.  In  order  to  most  effectively  develop  a  numerical 
approximation  the  Water  Resources  Research,  Journal  of  the  Hydraulics 
Division,  Journal  of  Waterways  Harbors  and  Coastal  Engineering  Division, 
International  Journal  for  Numerical  Methods  In  Engineering,  Advances  in 
Water  Resources,  and  Applied  Mathematical  Modeling  were  searched  over 
the  period  of  holdings  at  the  Waterways  Experiment  Station  Library. 
Additional  references  were  also  obtained  from  the  Journal  of  Computa¬ 
tional  Physics. 

The  numerical  approximation  of  the  transport  equation  is  an  active 
research  area  within  each  of  three  major  numerical  analysis  disciplines: 
finite  difference,  finite  element,  and  method  of  characteristics.  Re¬ 
view  was  limited  to  finite  difference  techniques.  The  following  eight 
major  methods  presented  in  Table  I  were  investigated  and  are  reported  in 
turn  in  detail.  Additional  techniques  found  in  the  literature  are 
briefly  outlined  along  with  practical  considerations  in  selecting  a 
numerical  method. 

1 .  Sheng  Work  at  Aeronautical  Research  Associates  of  Princeton  (ARAP) 
Sheng  [4]  considers  tlie  following  equation 


9c  ^  9uc  ^  9wc  _ 
0t  9x  9z 


(1.1) 
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Table  I 


Major  Numerical  Methods 


Method 

1.  Centered  space 

differences  with 
ARAP  smoothing 

2.  Flux-Corrected 

Transport 


3.  Holly-Preissmann 
Scheme 


4,  Method  of  Second 
Moments 


5.  Balanced  Expansion 
Technique 


6,  Stone  and  Brian 

Technique 

7.  Gray  and  Binder 


8.  Leendertse  Scheme 


Advantages 

Consistcul  With  present 
hydrodynamics  (ADI) 

Positivity  of 
concentrations 


Explicit  (coding  may 
be  simpler) 


Explicit  (coding  may 
be  simpler) 


Phase  error  may  be 
systematically 
reduced 


Consistent  with  present 
hydrodynamics  (ADI) 

Consistent  with  present 
hydrodynamics  (ADI) 

Consistent  with  present 
hydrodynamics  (ADI) 


Disadvantages 

Potential  for  negative 
concentrations 


Two  schemes  must  be 
implemented  (higher 
and  lower  order) 

Explicit  time  step 
limit.  Not  devel¬ 
oped  for  two 
dimensions 

Diffusion  not  consi¬ 
dered.  Explicit 
time  step  limit. 

Explicit  and  diffusion 
time  step  limit. 

Not  developed  for 
two  dimensions. 

Potential  for  negative 
concentrations 

Potential  for  negative 
concentrations 

Potential  for  negative 
concentrations 


The  following  finite  difference  form  of  (1.1)  is  considered. 


n+1  _  n  _  At  ,, 
'"ijk  ^i,k  AxAz  R 


(1 


where 

fj^  H  right  flux  for  grid  cell  surrounding  c^ 

f  =  left  flux  for  grid  cell  surrounding  c.  . 

1j  1  y  iC 

g  E  top  flux  for  grid  cell  surrounding  c.  , 

J.  X  ^  K 

gg  =  bottom  flux  for  grid  cell  surrounding  c^ 

The  quantities  in  (1.2)  are  indicated  in  Figure  1.  Several  approxima 
tions  are  presented  for  f^  ,  f.  ,  ♦  and  g^  based  upon  several 

alternative  difference  schemes.  These  are  presented  in  turn  below. 


UPWIND  SCHEME 


Let  us  consider  the  case  u.  ,>0,  u.  ^,>0  and  w.  , 

>0  to  illustrate  the  final  form  by  substituting  the  above  re¬ 
lations  under  these  conditions  into  (1.2)  above: 


n+1  n  .  . 

c.  ,  -  c.  ,  (u.  ,c,  ,  -  u.  ,  ,  c.  ,  ,) 

i,k  i>k  1 ,  k  1 ,  k _ 1-1  k  1-1. k 


(l.A) 


(w.  ,C.  ,  -w.  ,  ,C.  ,  t) 

i*ki*k  1, k“ 1  1 , k- 1 


If  u.  ,  =  u.  ^  ,  and  w.  =  w.  ,  ,  ,  then  the  above  equation  reduces 

to  the  following  form: 


n+1  n 
c.  ,  -  c.  , 
i,k  i,k 


^i.k  _  "^l.k  _  . 

Ax  ^i,k  ^i-l,k^  Az  ^^i,k  ^i,k-l 


(1.5) 


and  the  space  differences  are  only  first  order  accurate.  In  computation, 
an  artificial  (numerical)  dispersion  is  introduced. 

CENTRAL  DIFFERENCE  SCHEME 


r  ^^i+l,k  ^i,k^  , ,  . 

■  “l,k  - ^ 


r  *°l-l,k  'i,k' 

■  “l-l,k  -  2 - 


(1.6a) 


(1.6b) 


=  '^i.k  - - 


(1.6c) 


=  '"i.k-l 


^‘'i  k  "^i  k-1^ 


(1.6d) 


Substituting  the  above  relations  in  (1.2)  obtain 
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This  scheme  is 


Analogous  expressions  hold  for  g  and  g  . 
similar  to  second  upwind  differencing  and  uses  the  central  difference 
approximations  as  often  as  possible.  Lower  order  differencing  (first 
upwind)  is  employed  as  necessary  to  eliminate  negative  concentrations. 
The  effective  advective  concentration  employed  in  an  outgoing  flux  is 
always  less  than  or  equal  to  the  concentration  of  the  cell  providing  the 
flux. 

CENTRAL  DIFFERENCE  SCHEME  WITH  SPATIAL  SMOOTHING 

A  smoothing  scheme  developed  by  Sheng  was  employed  to  the  central 
difference  scheme  (1.7).  It  appears  in  most  instances  negative  concen¬ 
trations  may  be  prevented.  However,  it  is  highly  desirable  to  avoid  all 
negative  concentrations. 

The  smoothing  scheme  examines  the  solution  surface  in  both  coordi¬ 
nate  directions,  independently,  and  determines  whether  oscillations  are 
short  or  long  wave  phenomena.  Short  wave  oscillations  are  smoothed. 

The  smoothing  procedure  is  as  follows  for  v.  . 


where 


=  >1  - 


j|/ 

j  -"j-ll/ 


\  =  Iv.  - 


At  =  |Vj+T  -  ''j-i|/2Ax 
y  >  2 

then  the  following  curvature  check  is  performed. 
If 


2  2  2  2 
A  X  Aj,  <  0  or  A  x  <  0 
K  L 


(1.11) 


where 


V.  ,  -  2v.)/ax^ 

1-1  3  / 


i  -  * 

aJ  =  (v,^o  +  V,  -  2vj^^)yAx^ 


"R 


3+2 


\  ''j-2  - 


then  smoothing  is  applied  in  the  following  manner. 


V.  =  V.  +  6(v  +  V  -  2v.)  (1.12) 

3  3  3+1  3-1  3 

A 

where  is  the  smoothed  value  for  v^  and  3  is  a  positive  con¬ 

stant.  In  practice,  tests  have  suggested  y  =  4  and  3=1/4  are  the 
best  values  to  employ  in  this  smoothing  procedure. 
FLUX-CORRECTED-TRANSPORT  (FCT)  SCHEME 

This  method  was  originally  developed  by  Boris  and  Book  [6].  It 
has  been  subsequently  improved  and  generalized  by  Zalesak  [7].  It  is  a 
two  step  method,  first  involving  a  low  order  calculation  and  then  a  cor¬ 
rection  to  a  higher  order.  The  upwind  scheme  is  used  to  compute  the 
first  order  result 


39 


td  _  n  _  At 
i  ,  k  i ,  k  AxAz 


(tj  -  +  gj  -  g^  (1.13) 


where  c.  ,  is  the  first  order  (transported  and  diffused)  concentra- 
1 ,  k 

tion.  A  higher  order  scheme,  e.g.,  the  central  difference  scheme,  can 

2  2  2  2 

be  applied  to  compute  higher  order  fluxes  f^  »  fj  »  8-^  »  Sn  • 

K  1  •  1  D 

Antidif fusive  fluxes  are  then  defined  as 


^  =  4  -  4 
\  =  4  -  4 

A  -  2  1 


(1.14) 


^B 


It  is  these  antidif fusive  fluxes  which  are  limited  in  the 


Zalesak  [7]  procedure  such  that 


\  ’  °i+l/2,k 


0  1  W2,kl  ^ 


\  •  °i-l/2,k 

N'  "  ’  0i,k+l/2 


0  <  D.  ,  ,  <1 

-  1-1/2, k- 


0  1  1  ^ 


(1.15) 


\  \  '  Oi,k-l/2 


0l0i,k-l/2ll 


Finally 


n+1  _  td  _  At 
^i,k  ^i,k  AxAz 


(a^  -  -  a')  (1.16) 


The  determination  of  the  D  coefficients  will  be  presented  subsequently. 
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.  Fully  Multidimensional  Flux-Corrected  Transport  Algorithms  for  Fluids 


In  this  article,  Zalesak  [7]  develops  a  new  flux  limiter  which 
generalizes  to  multidimensional  problems.  The  new  flux  limiter  in  one 
dimension  is  shown  to  exhibit  superior  characteristics  over  the  original 


limiter . 


Zalesak  considers  the  two-dimensional  problem  in  the  following 


fashion.  Consider 


Ltd  n  ,  n'I/t-L  _L 

i.J  i,J  13  ij  V  1+1/2, j  1-1/2, j 


(2.1) 


+  F^  -  F^  ) 

i,j+l/2  '^i, 3-1/2/ 


Htd  n 

”i,j  ”i,j 


i-l/2,j 


(2.2) 


i,j+l/2  ■  ''i,j-l/2j 


where  ^  represents  the  lower  order  transported  and  diffused 

solution  and  w^*''f  the  higher  order  solution.  We  observe  that  the 
i  >  J 

difference  between  the  time  advancement  may  be  written  as 


Htd  Ltd  ..  .  ^-if/fH 

w.  .  -  w.  .  =  -  (Ax.  .Ay.  .)  (F. .1 /o  •  ~ 

i,:  i,3  1,3  1,3  L'  ^+1/2,3 


i+1/2, 
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Note  this  difference  is  written  as  an  array  of  fluxes  between  different 
grid  points  and  is  the  condition  required  to  implement  flux  corrected 
transport . 

Flux-corrected  transport  is  implemented  in  the  following  fashion 
using  implicit  difference  schemes  for  wV^f  and  w^^*?  : 

1.  Compute  from  (2.1)  above 

2.  Calculate  (w” 

\  i,J  i»J/ 

3.  Compute  wV^*?  from  (2.2)  above 

1,1 

/  r  1  1  r-H  /  n  Htd\ 

4.  Calculate  F  (w.  w.  .  ) 

\  1,1  1,1/ 

H  L 

5.  Determine  the  antidif fusive  fluxes  F  -  F  as  in  (2.3)  above 

6.  Limit  these  antidif fusive  fluxes  as  follows 


^i+l/2,j  ^i+l/2,j‘^i+l/2,j  ®  -  '^i+l/2,j  -  ^ 


^i,j+l/2  ■  ■^i,j+l/2^x,j+l/2  °  -  ^i,j+l/2  -  ^ 


(2.4) 


7.  Apply  the  limited  antidif fusive  fluxes: 

n+1  Ltd  .  .-1/  c  .c 

"l.J  ■  "l,j  -  (Vl/2,j  -  *1-1/2, J 

+  -  *li-l/2) 

The  crucial  step  in  the  above  process  is  Step  6,  the  flux- 
limiting  stage.  The  following  quantities  are  computed; 


P.  .  =  the  sum  of  all  antidif fusive  fluxes  into  grid  point  (i,j) 

1,1 


=  max(o,  \+l/2,p 


+  max(o,A.^._^/2>  '  ''i,j+l/2> 


(2.6) 
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""10+1/2 


(  ■  /„+  -  ^ 

''l,j+l/2  -  “  ) 

1,3+1’  ^1,3] 

V 

''l,i«/2  '  °  ) 

(2.13) 


Prior  to  determining  these  quantities  the  following  procedures  are 
employed . 


and  either 


/  Ltd  Ltd  \  .  r> 
(“i,j  -  ^  “ 


a  /  n  Ltd  \ 

w.  .  =  maxiw,  w.  .1 
iO  \iO  1,3/ 

b  .  /  n  Ltd\ 

w,  .  =  minlw.  . ,  w ,  .  1 
1,3  \1,3  1,3/ 

max  /a  a  a  a  a  \ 

w,  ,  =  maxiw.  ,  w.  w.,,  . ,  w ,  .  , ,  w ,  ] 

l,j  V  1-1,3  1,3  1+1,3  1,3-1  1,3+1/ 

min  .  /  b  b  b  b  b  \ 

i,j  V  1-1,3  1,3  1+1,3  1,3-1  1,3+iy 


(2.14) 


A  n  ^  A  /Ltd  Ltd\  ^  o 

Vl/2.j  ■  °  f”"  ''i+1/2.A“i«,J  ■  ”l.j)  ‘  ° 

and  either  j  "  "m,j)  “  ‘^.U) 

/  Ltd  td  \  ^  (s 
or  A. ..  .Iw.  .  -  w,  T  .  I  <  0 
1+1/2, j\  1,3  1-1,3/ 

A  n  C  A  /  Ltd  Ltd\  . 

and  either  -+2  ’  ““j+l)  ^ 


(2.16) 


(2.17) 


(2.18) 


(2.19) 


Zalesak  notes  that  while  the  solution  will'  be  kept  between  w 
and  w™^^  monotonicity  in  one  coordinate  direction  in  rare  cases  may  be 
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violated.  Additional,  more  stringent  limitations  are  given  to  prevent 
this  occurrence  but  will  not  be  presented  here. 

3.  Holly-Preissmann  Scheme 

In  developing  the  scheme  the  following  one-dimensional  advection 
equation  is  considered  for  constant  velocity  by  Holly  and  Preissmann  [8] 


(3.1) 


where 


c(x,t)  =  concentration  at  x  for  time  t 
u(x,t)  =  u^  a  constant  velocity 
For  this  case,  the  formal  solution  to  (3.1)  becomes 


c(x,t  +  t)  =  c(x  -  U^T,t) 


(3.2) 


Therefore  the  concentration  at  x  for  time  t  +  t  is  determined  from 
a  knowledge  of  the  concentration  for  time  t  at  x  -  u^t  .  This  forms 
the  basis  for  the  development  of  many  explicit  direct  calculation  methods. 
Holly  and  Preissmann  consider  the  situation  shown  in  Figure  2. 

In  order  to  determine  c^  an  interpolation  procedure  is  needed. 

In  this  scheme  information  at  only  the  two  adjacent  grid  points  is  used. 
Knowledge  of  c”  ,  employed,  where 

cx^  =  9c/3x  at  X  =  Xj^  ,  t  =  t^  . 

A  dimensionless  argument  a  ,  known  as  the  Courant  number,  is 
defined : 


^i  - 


(3.3) 


Letting  u  t  =  x  ,  a  general  distance  measured  from  x^  to  » 
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and 


(3.4) 


da  _  /  v-1 

^  -  ^i_i) 


dx  da  dx  da  i  *1-1"^ 


The  following  Interpolating  polynomial  is  considered. 


3  2 

y(a)  =  Aa  +  Ba  +  Da  +  E 


(3.5a) 


y(o)  =  c”  =  E  ,  y(l)  =  =  A  +  B  +  D  +  E  (3.5b) 


•,  ^  (3Aa  +  2Ba  +  D) 

,  -  v?  " 


(3.6a) 


i(0)  =  ex'"  _ _ 5 _  y(l)  =  (■14-+_^--t^).  =  cx" 

^  i  (x^  -  x^_^)  ’  (x^  -  x^_^)  i-1 


=  cx.  ^  (3.6b) 


Equations  (3.5b)  and  (3.6b)  are  employed  to  solve  for  A  ,  B  ,  D  , 
and  E  .  We  note  E  and  D  are  available  directly;  namely, 


E  =  c”  ,  D  =  (x^  -  x^_j^)  cx^ 


(3.7) 


From  y(l)  ,  we  obtain 


n  ,  .  n  n 

A  =  -B  +  c^_^  -  (x^  -  x^_^)  cx^  -  c. 


(3.8) 


Substituting  in  y(l)  ,  we  obtain 


«l_l(5Ci  -  x^_^)  =  -3B  +  3(c"_^  -  c")  -  3(x.  -  x^_^)  cx^ 


+  2B  +  (x^  -  cx^ 


(3.9) 


B  =  3(c;_^  -  cf)  -  (Zcxl  +  cx"_^)(x^  -  x._^)  (3.10) 
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Substituting  B  into  (3.8)  we  obtain 


A  =  -2(c”.^  -  c^)  +  (cx"  +  cxJ_J(x.  - 

Substituting  these  results  into  (3.5a),  we  obtain 
y(a)  =  |^(cx"  +  cx”_^^(xj^  -  x^_^) 

-  (2cx^  +  cx^.  J(x^ 

Collecting  terms  (3.12),  y(a)  is  rewritten  as 
y(a)  =  (-2a^  +  3a^) 


-  +  CX^(x^ 


-  X .  ,  )a  +  c . 
1-1  1 


(3.11) 


(3.12) 


+  (2a^  -  3a^  +1)  c"  +  (a^  -  a^) (x^  -  cx”_^  (3.13) 

+  (a^  -  2a^  +  a)(x^  -  cx” 

Thus  we  have 

=  y(a)  (3.14) 

In  order  to  advance  the  solution  in  time  the  concentration  derivative 
also  must  be  determined.  For  the  constant  velocity  case,  differentiating 
the  transport  equation  with  respect  to  x  and  interchanging  the  order 


of  the  X  and  t  derivatives 


^■."^■^.".  4'!  -V  v'\  ^yv  \  .".  i""yiiiiiii  ■  1 1 


K: 


‘1'- 


Ml 


l9 

* 


Using  cx  =  3c/3x  the  equation  above  may  be  written  as 


3cx  .  3cx  „ 

-r- - 1-  -r —  =  0 

3t  o  3x 


(3.16) 


We  note  that  this  form  is  exactly  analogous  to  that  in  (3.1)  for  the 
concentration  itself.  Considering  (3.6)  and  employing  previous  results 
for  the  A  ,  B  ,  and  D  coefficients 

y(a)  =  I  [3(x^  -  x._p(cx^_3^  +  cx")  -  6(c"_^  -  c")j 

+  -  c^)  -  2^2cx"  +  cx"_^)(x^  -  x^_j^)jot  (3.17) 


+  cx"  (X.  -  x._^)^  /(x^  -  x._^) 


Upon  rearrangement  we  obtain 


y(a)  =  3a^^cx"_j^  +  cx")  - 


(x.  -  X.  1 ) 
1  1-1 


+  itill - ^  a  _  zfzcTL^  +  cx"  a  +  cx" 

(Xi  -  x^_j^)  V  1  i-1/  1 


(3.18) 


or 


y(c)  -  /<■---  °\  c"  ,  +  c"  +  (3a^  -  2»)  cx" 

(x^  -  x^_j^)  i-l  (x^  -  x^_^)  1  i-1 


(3.19) 


+  (3a^  -  4a  +  1)  cx" 


Then  by  analogy 


n+1  • X  X 

cx^  =  y(a) 


(3.20) 
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A  linear  stability  analysis  is  performed  to  demonstrate  the  following 
stability  criteria 

^2  j  <  <1  for  a  £  1  (3.21) 

where 

jx^l  =  primary  mode 

1^21  -  ^  parasitic  computational  mode  sufficiently  small  relative 
to  unity  so  that  the  mode  disappears  from  the  solution 
quite  rapidly 

In  order  to  perform  the  computations,  initial  and  boundary  condi¬ 
tions  for  both  the  concentration  and  its  space  derivative  must  be  speci¬ 
fied.  The  specification  of  the  spatial  derivative  may  not  be  straight¬ 
forward  in  practical  computations. 

The  method  is  extended  to  the  non-uniform  velocity  case  in  the 
following  manner.  In  this  case  the  exact  solution  to  (3.1)  requires 
that  the  concentration  be  constant  on  the  trajectory  or  characteristic 
curve 


dx 

dt 


u(x,t) 


(3.22) 


Thus,  in  the  application  of  (3.14),  the  interpolation  argument  corre¬ 
sponds  to  the  point  on  the  x-axis  where  the  trajectory  defined  by  (3.22) 
crosses  it.  This  point  is  estimated  by  means  of  a  suitable  integration 
of  u(x,t)  from  point  i  to  i-1  over  the  time  interval  t  .  In 
order  to  compute  the  advection  of  the  concentration  derivative,  equa¬ 
tion  (3.16)  is  written  more  generally  as 


3cx  .  ,  3cx  ,  9u(x,t)  rt 


(3.23) 
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The  first  two  terms  are  evaluated  using  the  trajectory-estimated  a 
in  the  application  of  (3.20).  The  contribution  of  the  third  term  is 
determined  by  averaging  cx  and  3u/3x  from  the  point  i  to  i-1  . 

The  extension  of  the  technique  to  two  dimensions  for  Lhe  constant 
velocity  case  is  briefly  outlined.  However,  the  treatment  for  the  non- 
uniform  case  is  not  completely  apparent.  Diffusion  may  be  incorporated 
within  the  interpolation  scheme  but  only  at  the  expense  of  a  more  severe 
restriction  on  the  time  step  than  a  £  1  .  By  decoupling  the  diffusion 
calculation'^  from  advection  the  authors  suggest  this  restriction  may  be 
avoided. 

4 .  Method  of  Second  Moments 

Egan  and  Mahoney  introduced  this  scheme  in  meteorological  studies 
(air  pollution  transport)  [9].  Width  correction  adjustments  were  re¬ 
ported  by  Pedersen  and  Prahm  [10].  Pepper  and  Baker  [11]  have  developed 
an  elaborate  three-dimensional  transport  algorithm  for  predicting  tritium 
releases  from  the  Savannah  River  Nuclear  Power  Plant. 

The  basic  method  involves  describing  the  concentration  distribu¬ 
tion  within  each  cell  of  an  Euleiian  mesh  by  its  first  three  moments, 
zeroth  (total  mass),  first  (mass  center),  and  second  moment  (variance). 
The  cell  distribution  representations  are  then  advected  using  the  veloc¬ 
ity  components  developed  from  a  separate  solution  scheme.  At  the  end  of 
the  time  step  the  resulting  individual  distributions  are  combined  in  a 
composite  and  the  process  is  repeated  over  the  subsequent  time  steps. 

The  scheme  is  explicit  and  quasi-Lagrangian.  Due  to  the  explicit  nature 
of  the  advection  scheme  the  particle  Courant  number  must  be  less  than  or 
equal  to  one  to  maintain  stability. 
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The  method  is  presented  in  one-dimension  by  considering  a  concen¬ 
tration  distribution  to  be  represented  by  three  parameters  defined  as 
follows 

40.5 

C  =  f  C(£)  de 
m  J  m 

-0.5 

40.5 

C  F  =  f  C(e  )  e  de  (4.1) 

mmj  mm 

-0.5 


4-0.5 

C  =12  f  C(e  )(e  -  F  de 

mm  J  m  m  m 

-0.5 

where 

Ejn  =  relative  displacement  of  material  in  the  mth  cell  relative  to 
the  cell  center 

C  =  mean  (cell)  concentration  of  the  distribution 
m 

F  E  center  of  mass  of  concentration  distribution 
m 

2 

R  E  scaled  second  moment  of  concentration  distribution 
m 

2 

R  is  scaled,  multiplied  by  a  factor  of  12,  such  that  a  rectangular 
m 

2  2 

distribution  of  length  L  has  R  =  L  .  In  applying  the  method, 

m 

rectangular  type  distributions  are  maintained  in  each  grid  cell.  Ini¬ 
tially,  C  is  specified  and  the  distribution  is  assumed  to  occupy  the 
m 

entire  cell.  For  rectangular  type  distributions,  the  integrals  in 
(4.1)  above  are  evaluated  by  considering  for  each  grid  cell  the  material 
distribution  remaining  in  the  cell  and  those  which  entered  during  the 
time  interval. 

Figure  3  illustrates  the  procedure  for  advectlon  in  one  time  step 
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of  a  =  uAt/Ax  .  A  proportioning  parameter  is  defined  as  =  (F^  +  a 

+  R  /2  -  0.5)/R  .  For  P  <0  none  of  the  material  is  advected  into 

m  m  m  — 

cell  nri-1  .  For  P  >1  all  of  the  material  in  cell  m  is  advected 

m 

into  cell  nri-l  .  For  0<P  <1,  PC  is  advected  into  cell  m+1  , 

m  m  m  ’ 

while  (1  -  P  )C  remains  in  cell  m  .  Thus  in  general,  one  obtains 
mm  e>  » 


T+1 

C  =  C  +  C 
m  r  a 


T+1  T+1 

C*'**"?^'**"  “  C  F  ^  C  F 
mm  r  r  a  a 


(4. 2a 
(4.2b: 


cT+1  2  T+1 
m  m 


=  C  [r^  +  12('f^‘''^  -  F  +  C  [*R^  +  12 (f^'*'^  -  F  Vl 
rlr  \m  r/J  a|_a  \m  a/J 


(4.2c; 


where  subscripts  r  and  a  indicate  quantities  remaining  and  newly- 
advected  into  cell  m  ,  respectively. 


For 


P  <  0  ;  C  =  ,  F  =  F^  +  a  ,  R  =  R^ 

m  r  m  ’  r  m  ’  r  m 


P  >1:  C  =  0,  F  =  0,  R  =  0 
m  r  r  ’  r 


(4. 3a] 
(4.3b; 


C  =  (1  -  P  )C 
r  mm 


(i-r! 


+  P 


0  <  P  <  1  : 
m 


0  <  P  ,  <  1  : 
m-1 


C  =  P  ,C  , 
a  m-l  m-l 


(p  ,  - 
-  \  m-l  m-“l 


(4.3f) 


R  =  P  .R" 
a  m-l  m-l 


A  step-shaped  distribution  will  be  advected  downwind  without 
change  of  shape.  Small  diffusive  errors  will  remain  when  more  com¬ 
plicated  distributions  are  advected. 

The  technique  may  be  extended  to  two  dimensions  and,  in  the  fol¬ 
lowing,  we  employ  the  notation  of  Pedersen  and  Prahm  [10].  Two  por¬ 


tioning  parameters  are  defined  as  follows. 


‘^x>  -  ~-2  - 

r  ••  I 

Py  -  |sign  (cJyXFy  +  Oy) - ^y 


(4.4a) 


(4.4b) 


These  concepts  are  illustrated  in  Figure  4.  In  the  most  common  case, 

0  <  P  <  1  >  and  0  <  P  <  1  ,  and  the  computational  formulae  analogous 
X  y 

to  the  one-dimensional  case,  are  given  as  follows. 

Con t rlbutions  in  grid  cell  (m,n) 


<C™)^  .  <1 


P  )(1  -  P  )c 

X  y 


T  T  \ 

R  +  P  R  ) 

X  X  X  /  ^ 

- 2 - 


T  T 

R  +  P  R 


P  )r''’ 

X  X 

P  )r'^ 

y  y 


-  sign  (Oy) 


(4.5) 


55 


Figure  4.  Method  of  moments  advectlon 
procedure  (two  dimensions) 
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Contributions  in  cell  (uri-l  sign  (o  ) ,n+l  sign  (a  )) 

X  y 


PPG 
X  y 


(F 

X  a 


T+1 
(F  ^ 

y  a 


P  R  -  1) 

^  ""  f  -  sign  (o^) 


(Vy  -  0  - 


-  sign  (Oy) 


(4.8) 


T+1  T 

(R  )  =  P  R 

X  a  XX 


(R  =  P  r"^ 

y  a  y  y 


When  P^  and/or  P^  are  greater  than  1  and/or  less  than  zero,  differ¬ 
ent  equations  analogous  to  the  one-dimensional  case  must  be  applied.  We 
note  for  an  arbitrary  cell  (m,n),  3  distributions  may  be  advected  into 
the  cell  and  1  may  remain.  Formulae  analogous  to  (4.2)  are  employed  to 

compute  (C)  ,  (F  )  ,  (F  )  ,  (R  )„  „  »  and  (R^)_  „  • 

m,n  X  m,n  y  m,n  x  m,n  y  m,n 

Pedersen  and  Prahm  [10]  also  limit  the  width  of  the  distribution 
such  that  it  must  fall  within  one-cell  after  the  combination  process 
in  (4.2)  is  completed.  Analogous  limiting  may  be  performed  in  both 
coordinates  for  the  two-dimensional  case. 

This  technique  may  be  extended  using  fractional  steps  to  include 
the  diffusion  process  as  performed  by  Pepper  and  Baker  [11]. 

5.  Balanced  Expansion  Technique 


Chan  [12]  has  developed  a  new  procedure  to  construct  accurate 


finite  difference  advection  schemes  which  are  neutrally  stable 


(|x|  -l). 
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By  applying  the  procedure  In  a  systematic  manner,  the  phase  error  can 
be  reduced.  In  his  development,  the  following  model  equation  is 
considered. 

where 

({)  5  passive  scalar 
X  H  space  coordinate 
t  E  time 

u  >  0  and  is  a  constant  velocity 

A  Lagrangian  approach  is  employed  by  noting  that  the  value  of  (|)  is 
preserved  along  the  characteristic  dx/dt  =  u  in  (x,t)  space.  The 
general  solution  is  i|)(x,t)  =  F(x  -  ut)  ,  which  leads  to  the  following 
equations  in  the  discretization  of  (5.1). 


(|i(x,t  +  6t)  =  <t)(x  -  u6t,t) 
((i(x  +  u6t,t)  =  <ji(x,t  -  St) 


(5.2) 


where 

X  =  j6x 
t  =  ndt 

6t  =  time  Increment 
6x  =  space  Increment 

The  Courant  number  a  =  u6t/6x  is  introduced  and  the  procedure  is 
illustrated  in  Figure  5.  The  positions  of  the  six  quantities  > 

j  1  ’  1  ’  ^  symmetric  about  (|)^  ,  which  is 

midway  between  j-1  and  j  .  Each  of  these  six  quantities  is  expanded 
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in  a  Taylor  series  about  ({>^  .  The  following  balanced  differences  are 
constructed.  In  this  manner,  all  even  derivatives  in  the  Taylor's 
series  are  eliminated. 


■  I)  It  (“  + 1) 


^  4)  ^ 


<i)  +  It  (i)' 


ird)'  - 


“)  ''"♦x  *  It  (2  -  “) 


2/1  -  5^(5)  . 

IT  ( 2  -  ^  + 


From  (5.4)  let  us  solve  for  '5x<fi^  , 


(5.3) 


(5.4) 


(5.5) 


-  jr(j)  'hil) 


(5.6) 


Substituting  (5.6)  in  (5.3)  we  have 


♦"!i  -  (2<-  +  -  ♦"-!  -  It  (1/  ‘’■'♦xxx 

-  It  (i/  +  lr(“  1/  ‘-^♦xxx  +  lr(“  i) 


(5.7) 
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(1  +  2a) (1  +  a)  \  j-1 


,  (a  +  1)  (1  -  2a)  +  (2a 


(a  +  1) 


3  '^j-1 


(2a  - 
(2a  + 


t  ^j(a  ;  1)  [(2»  +  l)(l/  -  (»  +  l/] 


+  (1  -  2a)(|- 


2  X  5.5 

6x  (j) 


We  continue  to  simplify  by  noting 


(a  +  1)(1  -  2a)  +  (2a  -  1) (a  -  1)  =  -2a^  -  a  +  1  +  2a^  -  3a  +  1 


=  -4a  +  2  = 


and  the  6x^({)  term  may  be  written  as 


[at  -  (■  *  i)‘]  *  <‘  -  ->b)‘  -  a  -  •)*]  I  ^ 

■  ;.‘l »'  [©*  -  ("  *  t)1  *  '>  -  =•><“  ♦  -  (i  -  •)*]  j  ^ 


=  4.  q.  -  2a)  (1  -  g)  /_^n+l  _  ^n-1 


(1  +  2a) (1 


^  w-i  -  *r) 


2(1  -  2a)  /  n  n  \  ,  f 

+  1)  bi  ■  ^i-l)  + 


Dropping  the  truncation  error  a  discrete  formula  for  calculating 
is  developed  without  attempting  to  evaluate  derivatives.  The  computa¬ 
tional  scheme  in  (5.14)  is  a  three  time  level  scheme  containing  informa¬ 
tion  at  only  two  space  points  (j  and  j-1).  The  scheme  is  explicit  if 
computations  are  performed  sequentially  in  the  downstream  direction. 

Chan  obtains  several  other  formulae  by  using  balanced  pairs 
(centered  about  ^11  cases,  stability  is  governed  by  0  _<  a  ^  1  . 

By  including  more  points  symmetric  about  <|i^  ,  it  is  possible  to  develop 
formulae  which  successively  eliminate  higher  order  odd  space  deriva¬ 
tives.  In  (5.14)  4)^  and  have  been  eliminated.  Chan  presents 

schemes  which  eliminate  (j)^  as  well,  along  with  two  additional  formulae 

in  which  6  is  eliminated.  The  essential  characteristics  of  the 
XXX 

balanced  expansion  scheme  are:  (1)  all  even  space  derivatives  are 
eliminated,  which  Chan  notes  is  sufficient  to  insure  |x|  =1  ,  and 
(2)  successive  elimination  of  higher  order  space  derivatives  reduces  the 
phase  error. 

Chan  modifies  (5.14)  for  the  case  of  diffusion;  namely. 


ii  +  It  v 

3t  ^  3x  ^2 
3x 


where  v  is  a  diffusion  coefficient.  Figure  6  is  employed  in  develop¬ 


ing  the  method.  Unlike  in  Figure  5, 


is  not  equal  to 


which  is  the  value  of  <J)  located  at  a6x  to  the  left  of  <()j  .  He 


notes  that  the  fluid  particle  originally  at 


does  move  to  position 


at  the  end  of  the  6t  increment;  however,  the  diffusion  process 


has  changed  the  (j)  value  associated  with  the  particle.  Using  a  forward- 

2  2 

time,  central  space,  finite  difference  scheme  for  D<t)/Dt  =  v5  :)/3x  , 
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D/Dt  is  the  Lagrangian  derivative,  Chan  obtains 


(5.16) 


where  6  =  v6t/6x  .  For  stability  0  £  g  £  1/2  .  Similarly,  because 
of  diffusion,  tjij  ^  is  not  equal  to  •  They  are  related  by 


^  ^  (5.17) 


In  analogy,  with  the  original  development  for  advection  only,  the  six 

7n+l  ,n  Tn-l  Tn+l  n  ,  Tn-l 

quantities  ((»^_^  , 

panded  in  Taylor  series  about  cj)^  to  give 


(1  -  2a) (1  -  g)  /;n+l  'n-l\ 
(1  +  2a)  (1  +  a)  Tj-l  ~  ^3  / 


2(1  -  2a) 
(a  +  1) 


(5.18) 


The  solution  procedure  consists  then  of  the  following  three  steps 


1. 

^  7n-l 

Compute  (Pj 

using 

(5.17) 

for  the  entire  space  domain 

2. 

Use  (5.18)  to 

compute 

Tn+l 

♦j 

over  the  space  domain 

3. 

Use  (5.16)  to 

compute 

n+1 

♦j 

The  techniques  presented  and  preliminary  testing  by  Chan  show 
that  the  advection-dif fusion  schemes  are  very  accurate.  Unfortunately 
only  the  one  dimensional  case  is  considered  with  constant  coefficients. 
Extensions  to  multi-dimensional  problems,  flows  with  non-uniform 
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velocity  and  diffusion,  and  mesh  with  variable  spacings  have  not  been 
reported. 

6.  Stone  and  Brian  Technique  [13] 

In  their  development,  the  usual  one-dimensional  equation  is 
considered 

3x 

where 

D  =  diffusion  coefficient  (assumed  constant) 

V  E  constant  velocity 
u  =  concentration  of  a  given  material 
“  _  2  2 

Since  u(x,t)  =  T  A  e  sin  wif  (x  -  Vt)  satisfies  (6.1),  it 

w=l 

constitutes  a  solution  for  the  appropriate  boundary  conditions. 

A  general  form  of  finite  difference  approximations  to  (6.1)  is 
considered  as  follows 

-D  L^(u)  +  V  L^(u)  +  L^(u)  =0  (6.2) 

where 

2  _  2  2 
L^(u)  =  approximation  to  3  u/3x 

L^(u)  =  approximation  to  3u/3x 

L^(u)  =  approximation  to  3u/3t 

Stone  and  Brian  note,  that  the  corresponding  finite  difference  solution 
to  (6.2)  with  the  corresponding  boundary  conditions  may  be  written  as 
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~*.'  *.*  v’’  “7  -1*  •'* 


J-1 


u.  =  /  A  p’^  sin  WTf  (jAx  -  VAnAt) 
J  »n  JLi  w 

w=l 


(6.3) 


in  which  p  and  <()  depend  upon  wttAx  ,  VAt/Ax  ,  and  DAt/Ax  ,  and 

on  the  finite  difference  approximation  employed  and  J  represents  the 

number  segments  into  which  x  is  divided. 

It  is  their  objective  to  evaluate  the  accuracy  of  various  finite- 

difference  approximations  to  (6.1)  by  comparing  V  with  (jiV  and 
2  2 

-w  IT  DAt 

e  with  p  .  For  a  perfect  equation  (J)  =  1  for  all  J-1 

frequencies.  In  the  special  case,  D  =  0  ,  p  =  1  ;  however,  for 
D  2.  0  >  P  would  be  smaller  for  high-frequency  harmonies  than  for  low 
frequency  harmonies. 

The  general  finite  difference  analog  to  (6.1)  is  written  as 
follows. 

„  r.  2  ^“l.n  ■*'  "l.n+l'l  ,  V  r  ,  , 


+  ^  (u.  -  U.  1  )  +  c(u 

2  '  J ,n  j-l,n" 


j+l,n+l  ~  ^^^,0+1^ 


(6.4) 


‘^^''j,n+l  ''j-l,n+l^]  At  [^^^"3,0+1  "  *'j,n^ 


where 


a+'^+c  +  d  —  1 


69 


(6.7) 


a+-2  +  c  +  d=  l  =  2a  +  e 


g+-2+in=l  =  g  +  0 


Considering  the  above  relations,  and  D  =  0  ,  we  now  evaluate  (6.6) 


I  sin  wnAx 

BwirAx  I  (1  _  0)  +  +  0  _  00.  _  wttAx 


+  tan 


Y  sin  wttAx 


-  9)  +  s(i  -  9)  +  [e  -  e(|  - 


COS  wttAx 


2 

gwirAx 


sin  wnAx 


1  +  60  -  -  ej(l  -  cos  wttAx) 


(6.8a) 


(6.8b) 


If  one  takes  the  limit  of  (6.8b)  as  8  -»■  0  ,  noting  tan  "^x  :  x  ,  as 

X  ->■  0 


sin  (wttAx) 

lira  (fi _ (wttAx) _ 

B->-0  1-0  [1-  cos  (wttAx)] 


(wttAx)  ‘ 

-  0[-^j 


(wttAx)' 


(wttAx)  ^  (wttAx)  ^ 


(wttAx)  ^  (wttAx)  ^ 

A!  6! 


(6.9) 


If  0  =  1/3  ,  we  match  terms  through  the  second  order  and  ^  :  1 
Stone  and  Brian  note  that  e  =  1/2  is  the  best  value  for  the  single 
remaining  degree  of  freedom.  They  note  0  °  ^  ^  ~  ’ 
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0  =  1/3  .  Thus  6  =  1/3  insures  good  values  of  (j)  for  the  low  fre¬ 
quency  harmonies  (wtt  small)  as  6^0,  and  e  =  1/2  insures  that  these 
values  of  (j)  do  not  change  rapidly  for  a  range  of  B  slightly  greater 
than  zero. 

Collecting  all  previous  results  (e  =  1/2  ,  0  =  1/3  ,  a  =  1/4 

=  d  ,  c  =  1/4  ,  m  =  1/6  ,  g  =  2/3)  we  obtain  rewritting  (6.4)  the 
following  scheme. 


-D 


(u. -  u.  ) 
j+l,n  j,n 


'"j.n  -  <“3+1, n+1  '  "j.n+l’ 

^  4  <“j  ,n+l  “j-l,n+l^j  ^  at  ^3  <“j,n+l  “j 


(6.10) 


<“3-l.n+l 


u 


j-l.n 


j+l,n+l 


j+1 


.n^_ 


=  0 


The  interesting  feature  in  (6.10)  is  the  spatial  time  derivative  or 
"spread"  time  derivative.  If  we  also  return  to  (6.5),  we  now  obtain  for 
the  scheme  in  (6.10) 


2 

P 


cos  wir&x  -  a  sin^  ^w^Ax^  ^  (witAx)  ^ 

^  cos  wttAx  +  a  sin^  ^vmAx^j  ^  (wTTAx)j 


(6.11) 


for  all  wttAx  and  for  all  a  and  B  •  Therefore  the  method  is  uncon¬ 
ditionally  stable  with  no  time  step  restriction.  Stone  and  Brian  also 
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reason,  the  scheme  is  not  further  considered. 

Stone  and  Brian  consider  the  use  of  a  cyclic  set  of  difference 
equations.  If  N  different  finite-difference  equations  are  used  over 
N  time  steps  and  the  cycle  repeated,  then  in  (6.3) 


p  -  (pjPjj...p^) 


1/N 


<l>  = 


>1  +  +  ••• 
N 


(6.1A 


Stone  and  Brian  considered  N  =  3  and  developed  a  different  6  for 
e  =  1/2  in  each  of  three  finite  difference  schemes.  In  each  scheme 
p  =  1  for  D  =  0  and  the  previous  relationships  among  the  coefficients 
are  sufficient  to  completely  define  the  scheme  once  0  is  determined. 
Consider  (6.9)  as  follows: 


lim 

P  -»■  0 


lim 


(•tj) 


lim 

e  -»■ 


(tii) 


lim 


^  0 


(^III> 


(6.15) 


Letting 


I 

i 

I 


Equating  the  first  three  powers  in  (wit Ax)  three  simultaneous  equa¬ 
tions  are  obtained  in  terms  of  0^^  ,  0^  ,  and  0^  .  This  cyclic  use  of 

three  different  difference  equations  is  shown  to  be  superior  to  the 
multiple  application  of  the  single  equation  considered  previously. 

7 .  An  Analysis  of  the  Numerical 

Solution  of  the  Transport  Equation 

Gray  and  Finder  [14]  consider  the  one-dimensional  transport  equa¬ 
tion  with  constant  velocity,  written  in  their  notation  as 


3c  .  3c  ^  3^c 

dX 


(7.1) 


where 

c  =  constituent  concentration 
u  i  transport  velocity 
X  =  space  coordinate 
t  E  time  coordinate 
D  =  diffusion  coefficient  (constant) 

As  in  previous  work,  the  general  solution  to  (7.1)  is  considered  as  a 
Fourier  series. 

00 

c  =  2  c^  exp  +  io^x)  ,  where  |x  -  ut(  £  1  (7.2) 

n=-oo 

and  is  the  (time)  frequency  of  the  nth  component,  is  the 

spatial  frequency,  and  i  =  .  If  one  considers  a  single  component 

in  (7.2)  and  substitutes  it  into  (7.1),  the  following  relationship  is 


obtained 


3  +  ua  -  iDo  =  0 

n  n  n 


3  =  a  (iDo  -  u) 

n  n  n 


Thus  for  a  single  component  solution,  one  obtains 


c  ^  exp  [ia^(x  -  ut)]  exp  ^-Do^t^ 


(7.3; 


(7.4: 


(7.5) 


where  the  first  exponential  describes  the  translation  and  the  second 
describes  the  amplitude  modulation  of  a  Fourier  component  with  time. 
From  (7.2)  for  a  single  component 


=  c^  exp  [iBjj(t  +  At)]  exp  (io^x) 

=  c  exp  (13  t  +  io  x)  exp  (i3  At)  =  c^  exp  (13  At) 
n  ‘^  n  n  *^  n  t  ^  n' 


(7.6) 


Thus  exp  (13  At)  is  considered  an  eigenvalue,  X 
n  n 

Let  the  eigenvalue  of  the  nth  Fourier  component  obtained  from  the 
numerical  scheme  be  denoted  by  X^  ,  Considering  the  computed  and  ana¬ 
lytical  components  after  a  time  such  that  the  analytical  wave  has  propa' 
gated  one  wavelength,  the  number  of  time  steps  required  is  given  as 
follows . 


N 


n 


n  _  _n  Ax 
uAt  Ax  uAt 


(7.7) 


where  is  the  wavelength  of  the  nth  component.  The  ratio  of  the 

computed  to  actual  amplitude  after  one  wavelength  is  gi.*en  by 
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we  obtain  the  following  compact  relation. 


^  yx 
\3  3 


,  uAt  6x  DAt  .  2\  k+1 


(7.11 


Ax 


e)6xnc^ 


.  io  Ax 

If  we  consider  ~  ^i+1  k  ~  ^  Fourier  com¬ 

ponent  and  substitute  into  (7.11)  and  note  for  0  = 


.  2  „  1  -  cos  26  +  e-^2e  _  ^ 

Sin  0  =  - 2 -  "  - 14 


i0  -i0  i0  .  -i0 

6  6  616 

isin  0  =  - 2 -  ®  “  2 


we  may  write  directly  the  relation  for  . 


cos  o  Ax 

k+l  4  + - - (1  -  e)  [vl  sin  a  Ax  +  2D(1  -  cos  a  Ax)] 

,  ^  c _  ^  3 _ 3 _ _ n _ “ 

n  ”  k  "  -  cos  o  Ax  , 

^  -|  + - j-S —  +  c(vl  sin  o^Ax  +  2D(1  -  cos  o^Ax)] 


(7.12) 


For  stability  £  1  »  which  requires  1  ^  e  ^  0.5  . 

In  their  paper.  Gray  and  Finder  consider  the  following  finite  dif 
ference  approximation  to  (7.1)  written  in  compressed  operator  notation. 


uAt  6x  DAt  „  2\  k+l 

'  1  -  7T  )  ■= 

Ax  / 


-  fl  <1  -  c)  (1  -  O  6x^lc’' 


(7.13 


The  essential  difference  between  (7.11)  and  (7.13)  is  embodied 
in  the  treatment  of  3c/3t  .  In  (7.11)  a  spatial  weighted  "spread" 
time  derivative  approach  is  employed,  while  in  (7.13)  the  time  deriva¬ 
tive  is  not  spatially  weighted. 

The  corresponding  eigenvalue  for  (7.13)  is  written  below  as 
follows 

k+1  1  -  (1  -  E)[vi  sin  a  Ax  +  2D'(1  -  cos  a  Ax)3 

T.  '  =  ~ -  =  _ S _ H _  (-7 

n  1  +  e[vi  sin  o^Ax  +  2D’ (1  -  cos  o^Ax)]  '■  '  ' 

Again  |a^|  £  1  ,  for  0.5  _<  e  _<  1  . 

Gray  and  Finder  employ  equation  (7. 7-7. 9)  to  study  the  character¬ 
istics  of  the  two  schemes.  An  eigenvalue  amplitude  plot  for  e  =  0.5 
and  1  is  developed  versus  “  (2tt/L^)Ax  =  (27r/nAx)Ax  =  2iT/n  for 

n  £  2  for  both  (7.12)  and  (7.14). 

Equation  (7.8)  is  employed  to  provide  amplitude  ratio  plots  for 
both  schemes  for  e  =  0.5  and  1  ,  while,  Equation  (7.9)  is  used  to 
develop  phase  angle  plots  for  e  =  0.5  and  1  .  Based  upon  these  plots 
and  direct  numerical  simulation  of  a  step  function  concentration  dis¬ 
tribution  the  spread  time  derivative  scheme  is  shown  superior  to  the 
standard  time  derivative  formulation. 

8 .  The  Leendertse  Formulation 

Leendertse  [15]  considered  the  space  staggered  grid  illustrated 
in  Figure  7  in  applying  the  finite  difference  approximations  to  the 
following  equation  sets 
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3(HP)  .  3 (HuP)  3(HvP)  3_  /  3P\  \  y  3y  / 

at  3x  3y  3x  V  X  3x  /  3y 


-  HS  =  0  (8,1) 

A 


2  2  1/2 

3u  .  3u  .  3u  ,  .  3n  .  u(u  +  V  )  1  _s 


(8.2) 


2  2  1/2 

3v  ,  3v  ,  3v  ,  3n  ,  v(u  +  V  )  __  1  s 


(8.3) 


In  ,  3 (Hu)  3(Hv)  ^  Q 
3t  3x  3y 


(8.4) 


where 


t  =  time 

x,y  =  Cartesian  coordinates 
H  =  water  depth 

u,v  H  depth  integrated  velocities  in  the  x  and  y  directions, 

respectively 

f  2  Coriolis  parameter 

c  =  Chezy  coefficient 

p  =  density  of  water 
s  s 

T  ,T  =  surface  stresses  in  the  x  and  y  directions, 
x’  y 

respectively 

P  =  pollutant  concentration 
n  H  water  surface  elevation 


g  =  acceleration  of  gravity 
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I  '  >  "If « 


& 


D  ,D  =  dispersion  coefficients  in  the  x  and  y  directions, 

X  y 

respectively 

S  =  source  strength  (concentration/unit/ time) 

A 

We  note  in  the  above  equation  set  that  P  is  considered  an  arbitrary 
pollutant. 

The  finite  difference  scheme  is  ADI  with  the  time  step  At  being 
split  in  half  to  advance  the  solution  to  t  +  At  .  In  following  the  fi- 


i 

nite  difference  formulations, 

for  details  of  notation. 

In 

' 

lowing  finite  differences 

are 

P* 

x-sweep  t  ->■  t  +  At/2 

k*’ 

JAx  -»■  X 

lii 

kAy  y 

1 

• - 

nAt  -»■  t 

For  continuity  we  write 

2  /  n+1/2  n  \  ^  f/'v 


+  n 


n  ) 

1  n+1/2  / 

'j.kj 

•“j+l/2,k 

,n  \ 

n+1/2  I  1 

'j.V 

j-l/2,kj  2Ax 

n  ) 

j.kj 

'''j,k+l/2  ■ 

»  i-  1 

j+l/2,k+l/2  "j+l/2,k-l/2  ^j+l,k 


,k+l/2  *’j-l/2,k-l/2  '^j-l,k 


.k+1 


,1  ^ 
,k+l/2  j-l/2,k+l/2  '^j,k 


'’j.k-l)  ''j,k-l/2]  2Ay  ° 
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(8.5) 


i 


which  may  be  written  in  general  as: 


n 


-r  4-  T-  n+1/2  _ 

j-1/2  j-l/2,k  +  '’j.k  ^  ''j+l/2“j+l/2,k  ■ 


(8.6) 


For  u  momentum: 


/  n-l/2  n-1/2  \ 

—  K  fv  +  \  J-«-3/2,k  ~  “1-1/2.  J  n+1/2 

At  ^“j+l/2.k  "j+l/2.k;  ^“j+l/2.k 

[Centered  Difference] 

/  n-1/2  _  n-1/2  )  ^  .  /„n+l/2  n+1/2  n-1/2  n-1/2  \ 

V“j+l/2.k+l  “j+l/2.k-l,/  2Ay  ^  (,^J+l.k  -  ’'j.k  +  ''j+l.k  -  ''j.k  ) 


2Ax 


[Centered  Difference] 


[Tiae  averaging  at  n+1/2] 
Derivative  centered  at 
time  level  n 


4  /  n+1/2  .n-1/2  nF/  n-1/2  Y 

+  2  V  J+l/2.k  Vl/2.kjK“l+l/2.kj 


+  V 


.ry  .  -X,  ,-x,2 
(h-'  +  n  )(c  ) 


1/2 


-/•cy  j.  “Xn  X 
pCh-'  +  n  ) 


T®  -  0 


(8.7) 


where 


1  /n 

4  Vj+1, 


k+1/2  ''j,k+l/2  ''j,k-l/2  ^  Vl,k-l/2) 


—y  1 

^  "  2  ^^j+1/2, k+1/2  ^3+1/2, k-1/2^ 


i  (vx,.  - 

(Vi,k  +  ""j.k) 


-X  1  /  n 
c  =2  ‘ 


The  above  equation  may  be  written  in  the  following  form: 


n+1/2  1  n+1/2  ^  n+1/2  _  „n 

j  j,k  ^  ''j+l/2''3+l/2.k  "'j+l'^j+l.k  ■  ®j+l/2,k 


(8.8) 


For  the  constituent  of  concern  we  obtain:  (Note  i  has  been  dropped, 
since  we  are  considering  only  a  single  constituent.) 
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n+y2  n+1/2  _  pn  n 
1 ,k  1 ,k  1 ,k  1 ,k 


r  n  n+1/2  /„n+l/2  n+l/2\ 

[®''j-l/2,k'‘j-l/2,k\  j-l/2,k  ^j,k  ) 


(‘'^k-l  +  ‘'Lk)  -  “yLk+l/2''5.k+l/2(’’J.k+l  +  '■j.k)] 

,  1  ^  n+1/2  /pn+1/2  n+1/2 \  _  n+1/2 


2  (Ax) 


3  n+1/2 
''j+1/2 


/  n+1/2  _  pn+l/2\1  ^  _1 _  T  n  ^  n 

j+l,k  j,k  )\  2(Ay)^  L  ^j.k-] 


«  (-5.k  -  "^k-l)  -  -  ^S.k)]  =  0 


where 


®’'j+l/2,k  (’^j+l,k  '’j.k  ^j+l/2,k-l/2  '^j+l/2,k+l/2) 

®yj,k+l/2  "  (’^j,k+l  '^j.k  ''j-l/2,k+l/2  ’^j+l/2,k+l/2) 

This  equation  is  then  written  ir  the  following  general  form: 


(8.9) 


a  +  b.p"+J/^  +  c.P"^y^  =  D. 

j  J-l,k  j  j,k  j  j+l,k  2 


(8.10) 


Let  us  now  consider  the  y  sweep  in  which  the  solution  is  advanced 
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from  time  level  n+1/2  ^  n  or  t  +  At/2  t  +  At  .  Leendertse  ex. 
presses  the  continuity  equation  as  follows: 


^+l/2t 


n  + 


+  A^)u]  + 


=  0 


j,k, n+1/2  (8.11) 


where : 


'^+l/2t'^ 


n+1  n 
'^j,k  ~  ^i,k 
At 


7h  4.  +  n+1/2  n+1/2  \  ' 

«  [(h>'  +  '^3-1/2, k  ^ 

X  X  2  2  / 


(h  4-  K  4.  n+1/2 

_v  j+l/2,k+l/2  '’j+l/2.k-l/2  '’j+l.k 


^  n+l/2\  n+1/2  /, 

"^j.k  )  '*j+l/2,k  ■  V^j-1/2, 


k+1/2 


.  .  ,  n+1/2  ^  n+1/2 \  n+1/2  ]  A. 

^j-l/2,k-l/2  '^j,k  "^j-l.kj  '^j-l/2.kj/^^’' 

- 

r  -I  /h  4-  h  n+1/2  \ 

[(h*  +  ^>')vj  .  6^  ^.j.+l/2,k  ^  1-1/2, k  ^  ^  "l.k-l/lj  ^ 


h. h.  ^  _n+l/2 


j+1/2, k+1/2  ^j-1/2, k+1/2  '’j.k+l 


4.  n+1  / 

^j.k  /  k+1/2  ■  (*^j+l/2,k- 


1/2 


+  h. 


j-l/2,k-l/2  '’j,k 


,  „n+l/2  , 

+  n,  ,  +  n 


n+1/2  \  n+1  1  L 

j,k-l/''j,k-l/2j/^^y 
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The  V  momentum  equation  for  the  y  sweep  is  given  as : 


gv 

(u)*-  +  (v_)‘-J 

-  J 

f  4.  ‘‘  \  r 

^.i.k+l/2  ^i,k+l/2)  ,=,2  ^  / 

(h^  +  n^) (c^)^ 

-  g  \ 

2  — /  y 

(N+l/2,k+l/2  ^i-l/2.k.+l/2 

2 


+  n 


n+l/2 

j.k+1 


+ 


n+1/ 


'j.k 


X 


■(: 


+  c  . 


s 

T 

_P 


s 

T 

-JL 

__E_ 


,^+l/2,k+l/2 


^.i-l/2>k+l/2 


n+l/2 

j.k+1 


+  n 


n+1 


The  constituent  equation  for  the  y  sweep  is  expressed  as: 


6+t/2tP(h  +  n)]  +  6^[(hy  +  n'')up"'']  +  6y[(h''  +  n^)v^py] 

-6^,[(hy  +  n'')D^6^pJ  -  6y[^(h^  +  ny)Dy  6yPj.J  =  0  at  j,k,n+l/2 


(8 


where 


'5+t/2tP(h  +  n)] 


-  C(^  ^ 


n+l\ 

j.k/ 


At 

2 


with 


~  ‘’j+l/2,k+l/2  ^j+l/2,k-l/2  ^j-l/2,k+l/2 

+  h 


‘j-l/2,k-l/2 


6  [(h^  +  n^)uP^]  =  6 

X  X 


,k+l/2  N.k-1/2 


^  n+1/2  ^  n+1/2  \ 

^+l/2.k  ^-l/2.kj 


'pn+1/2 

Zl+1/2 


+  \ 


‘j+l/2,k+l/2 


.  u  .  n+1/2  ^  n+l/2\  n+1/2  /„n+l/2 

^j+l/2,k-l/2  '^j+l.k  "^j.k  )“j+l/2,kVj+l,k 


■jtk' V"''*]  ■  [("3-1/2. 


k+1/2  ^j-l/2,k-l/2 
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(*'l-H/2. 
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r  -\  ( h  ,  +  h  4-  X  n+1/2  N 
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Note  that  the  above  equation  may  be  written  as 
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a,  .  P .  ,  ,  +  b,  P ,  ,  +  c,  ,  P .  -  =  D, 

k-1  j,k-l  k  j,k  k+1  j,k+l  k 

In  this  approach  the  constituent  equation  is  solved  directly  within  the 
hydrodynamic  computation  sweeps  although  there  is  no  coupling. 

In  the  application  to  Jamaica  Bay,  salient  simulation  character¬ 
istics  are  as  follows: 

a.  As  =  Ax  =  Ay  =  500  ft  (15.24  m) 

At  =  120  sec 

b.  The  explicit  time  step,  At^  ,  is  given  as  follows. 

4t  -  — — —  -  -  12.5  sec 

max 


c.  The  gravity  wave  Courant  number  Cr^  is  given  by 


Or  =  ^  =  10 
w  At 

e 


d.  The  particle  Courant  number  Cr^  is  given  by 


Cr  -  ^  .  (°.=  5  .  (,  ,,^20 

p  's  500  ft** 


e.  The  dispersion  coefficient  formulation  is  given  by 


D  =  14.3/2g  uHc  ^  +  D  ,  where  D  e(25,  45)  ft^/sec^ 

X  w  w 


where 


=  dispersion  coefficient 
g  =  gravity 


*  0.1524  m/sec. 

**  152.4  m. 

t  232,4.18  m  /sec. 
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V 


c  =  Chezy  coefficient 
H  =  local  water  depth 


£  background  dispersion  coefficient 
The  Chezy  coefficient  is  given  by 


where 


1.49  ,  1.49  . 

■  ^^1/6  '  (0.01)  (11) 


c  E  Chezy  coefficient 
R  =  hydraulic  radius 
n  E  Manning's  roughness 
g.  The  dispersion  coefficient  used  becomes 


0  =  +  25  =  37.5  ft^/sec  [3.484  m^/sec] 

X  100 


^  (37.5)  (120)  ^ 

.  2  250,000  ^ 

As 


and  a  cell  Peclet  number,  Pe  =  ~ 

X 

9.  Additional  Methods  and  Considerations 

Bram  van  Leer  [16,  17]  has  developed  several  upstream  centered 
higher  order  convective  schemes.  His  work  is  highly  theoretical  within 
the  domain  of  numerical  analysis  but  seems  to  indicate  a  general  ap¬ 
proach  to  constructing  extremely  accurate  finite  difference  schemes. 
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Forester  [18]  has  presented  a  non-linear  filtering  technique  for  higher 
order  even  (greater  than  two)  finite  difference  schemes  which  preserves 
the  peak  of  an  external  type  distribution  unlike  flux  corrected  trans¬ 
port  filtering.  The  application  of  the  filter  to  points  near  the 
boundary  was  not  reported  and  the  determination  of  the  coefficients  must 
be  made  through  direct  numerical  experimentation  for  each  problem 
individually. 

Narayan  and  Shankar  [19]  have  employed  a  multi-'=operational  scheme 
similar  to  the  Leendertse  scheme  previously  outlined  in  an  application 
to  Galveston  Bay.  Oster  et  al.  [20]  employed  Leendertse' s  scheme  [15] 
with  upwind  differencing  in  a  two-dimensional  computation.  Hinstrup 
et  al.  [21]  at  Danish  Hydraulic  Institute  have  developed  a  two- 
dimensional  explicit  scheme  employing  12  point  Everett  interpolation. 

The  treatment  of  boundary  cells  appears  to  result  in  some  mass 
falsification. 

Runchal  [22],  Siemieniuch  and  Gladwell  [23],  and  Lillington  and 
Shepherd  [24]  demonstrate  the  oscillatory  nature  of  central  difference 
approximation  to  the  steady  state  eq  .ation  in  convection  dominated  prob¬ 
lems  for  cell  Peclet  numbers  greater  than  2.  Jensen  and  Finlayson  [25] 
note  this  oscillatory  behavior  may  be  observed  in  transient  simulations 
of  sharp  fronts  as  well  for  improper  time  and  space  scales. 

Molenkamp  [26]  has  provided  a  review  of  several  finite  difference 
approximations  and  compared  them  for  a  rotation  of  a  circular  distribu¬ 
tion.  This  test  problem  provides  for  a  non-uniform  velocity  field 
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within  two  dimensions  and  will  be  employed  as  a  test  case  for  the  to  be 
developed  salinity  algorithms. 

The  cell  Peclet  number  as  mentioned  above  plays  a  major  role  in 
characterizing  the  relationship  between  the  grid  resolution  and  the 
numerical  scheme.  It  is  defined  as 


Pe  = 

X 


Pe  =  ^ 

y  K 

y 


(9.1) 


where 

Ax  =  grid  spacing 

u  =  maximum  magnitude  of  the  velocity  in  the  x  direction 
V  =  maximum  magnitude  of  the  velocity  in  the  y  direction 
=  dispersion  coefficient  in  the  x  direction 
Ky  H  dispersion  coefficient  in  the  y  direction 
The  Peclet  number  limit  of  two  is  required  to  prevent  oscillations  in 
the  solution  in  the  vicinity  of  a  sharp  concentration  front  for  central 
space  differencing.  For  typical  velocities  and  dispersion  coefficients. 
Ax  and  Ay  would  be  in  the  scale  of  hundreds  of  feet.  This  space 
scale  is  too  severe  to  be  applied  over  the  entire  area  of  Mississippi 
Sound.  The  Peclet  number  limit,  however,  is  only  significant  for  sharp 
fronts  and  although  there  may  be  significant  horizontal  gradients  in 
salinity,  these  gradients  are  not  as  severe  as  a  shock  or  discontinuity 
in  the  distribution.  If  first  order  upstream  differences  are  applied 
no  oscillations  will  develop,  but  accuracy  limitations  (such  as  those 
developed  by  Leonard  [27])  usually  require  a  dense  grid.  In  practical 
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computations,  one  normal  decides  on  a  space  scale  of  significance,  S  , 
and  selects  the  grid  spacings,  e  ,  such  that  e  _<  S/10  .  Higher  order 
schemes  may  allow  this  limitation  to  be  relaxed.  However,  these  higher 
order  methods  normally  involve  more  complicated  algorithms  and  increased 
computational  cost  and  model  development  time. 

In  conclusion,  it  should  be  noted  that  there  is  no  one  best  com¬ 
putational  finite  difference  scheme  for  the  transport  equation.  How¬ 
ever,  the  necessity  to  perform  computations  over  a  two-dimensional  grid 
with  irregular  boundaries,  suggests  that  a  simpler  lower  order  method 
be  selected  which  is  not  too  inaccurate. 


PART  IV:  NUMERICAL  METHOD  SELECTION  AND  DEVELOPMENT 


The  method  of  solution  for  the  transport  equation  must  be  compat¬ 
ible  with  the  hydr'^dynamic  scheme  employed  in  WIFM.  Since  the  convec¬ 
tive  terms  must  be  treated  in  the  hydrodynamics  for  tidal  circulation, 
the  following  explicit  type  time  step  limitation  must  be  obeyed 


uAt  <  min 


=  min 


/  Area  \ 
\ Diagonal / 


(3.1 


where 

u  =  maximum  particle  velocity  magnitude 
Ax  =  spatial  increment  in  x  direction  (variable) 

Ay  i  spatial  Increment  in  y  direction  (variable) 

It  is  desirable  to  leave  open  the  option  of  coupling  the  salinity 
transport  through  an  equation  of  state  to  the  density  involved  in  the 
pressure  gradient  terms  within  the  hydrodynamics.  If  this  is  to  be 
accomplished  the  method  of  solution  of  the  transport  equation  must  also 
satisfy  the  above  equation.  This  allows  both  explicit  and  implicit 
methods  for  solution  to  the  transport  equation  to  be  considered.  Ex¬ 
plicit  methods  must  obey  (3.1),  whereas  for  implicit  methods  the  only 
time  step  limit  is  one  of  accuracy. 

If  density  coupling  is  not  necessary  and  explicit  methods  are 

employed  then  it  may  be  possible  to  employ  a  time  step  in  the  transport 

T  H  H 

solution.  At  =  n  At  ,  where  n  is  an  integer  and  At  is  the  time 

step  in  the  hydrodynamic  solution.  In  this  case  the  following  limit 


must  be  obeyed  for  an  explicit  method 


-..T  .  .  /  AxAy 

uAt  <  min  I - ^ 


,  2  2 
Ax  +  Ay 


(3.2) 


where  u  represents  the  maximum  magnitude  of  the  particle  velocity  aver 

U 

aged  over  nAt  .  The  explicit  time  step  transport  limit  will  be  less 
than  an  implicit  time  step  transport  limit.  Thus  an  implicit  method  may 
accommodate  larger  time  steps  than  an  explicit  method.  Leendertse  [15] 
notes  the  problem  of  conservation  in  averaging  velocity  fields.  How¬ 
ever,  in  order  to  develop  long  term  transport  patterns,  it  may  be  desir¬ 
able  to  employ  much  larger  time  steps  in  the  transport  equation  solution 
than  in  the  hydrodynamics.  For  this  reason  and  to  maintain  consistency 
of  approach  with  the  hydrodynamics  a  multi-operational  implicit  scheme 
will  be  developed  for  the  transport  equation. 

As  has  been  previously  presented,  the  work  of  Stone  and  Brian  [13] 
and  Gray  and  Finder  [14]  illustrate  the  improved  computational  character 
Istics  of  the  spread  time  derivative  method  over  the  standard  forward 
time  centered  space  method  in  one-dimensional  problems.  Siemon  [28] 
has  investigated  the  extension  of  the  method  to  semi  two  dimensional 
problems  (advection  in  one  coordinate  direction,  diffusion  in  both 
coordinate  directions)  and  reported  favorable  results.  To  the  author's 
knowledge,  the  extension  to  completely  two  dimensional  problems  has  not 
been  reported.  It  is  proposed  that  this  extension  be  Investigated  in 
this  project.  By  employing  appropriate  coefficients  in  the  numerical 
formulation,  the  method  could  be  degenerated  to  Leendertse 's  approach. 


In  this  manner,  two  computational  schemes  may  be  coded  in  one  operation. 
The  format  of  the  spread  time  derivative  scheme  is  such  that  it  may  not 
be  expressed  in  a  form  suitable  to  flux  corrected  transport.  As  a  re¬ 
sult,  oscillation''  in  the  solution  may  be  smoothed  using  filtering 
techniques.  The  necessity  of  filtering  cannot  be  determined  until  nu¬ 
merical  experiments  are  conducted. 

As  an  alternative,  Leendertse's  approach  will  be  implemented  with 
flux  correction.  The  higher  order  scheme  will  correspond  to  the  stan¬ 
dard  Leendertse  formulation.  The  lower  order  scheme  will  employ  upwind 
differencing  of  the  advective  terms.  Thus,  two  schemes  must  effectively 
be  programmed  (a  higher  and  lower  order  scheme)  for  flux  correction. 

The  two  alternative  schemes  spread  time  derivative  and  Leendertse 
flux  corrected  will  be  compared  through  numerical  experimentation  to 
determine  the  most  appropriate  technique  for  application  in  Mississippi 
Sound. 

The  flux  correction  method  is  such  that  any  higher  order  method 

2  2  2 

may  be  employed.  Leendertse's  method  is  0(At  ,  Ax  ,  Ay  )  and  is  amenable 
to  adaptation  to  variable  grid  spacing  using  the  exponential  stretch 
transformation  in  WIFM.  In  the  future,  it  may  be  desirable  to  consider 
higher  order  compact  differencing  schemes  such  as  the  Kreiss  scheme  as 
reported  by  Roache  [29].  The  flux  correction  method  will  accommodate 
further  research  in  the  development  and  implementation  of  higher  order 
schemes. 

The  general  strategy  and  development  of  the  numerical  approxima¬ 
tions  to  the  transport  equation  for  application  to  Mississippi  Sound  is 


shown  in  Figure  8.  This  approach  provides  for  development  of  an  optimal 
second  order  method. 
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